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1.  INTRODUCTION 

\  i 

The  study  of  many  equilibrium  phenomena  leads  to  non-linear  equations 
which  involve  a  number  of  intrinsic  parameters.  Interest  then  centers  rarely 
on  the  determination  of  a  few  specific  solutions  of  the  equations  for  fixed 
parameter  values  but  rather  on  an  assessment  of  the  behavior  of  these  solutions 
under  general  variations  of  the  parameters.  For  example,  in  structural  analysis 
the  parameters  may  characterize  load  points  and  load  directions,  material  pro¬ 
perties,  geometrical  data,  etc.  The  set  of  all  solutions  and  associated  para¬ 
meter  values  has  been  called  the  equilibrium  surface  of  the  structure  (see 
eg.  [311).  This  equilibrium  surface  provides  considerable  insight  into  the 
behavior  of  the  structure  .and  the  stability  properties  (see  eg.  [23],  [32l 
for  further  discussions  and  various  examples).  From  a  numerical  viewpoint  the 
question  then  is  to  analyze  computationally  the  shape  and  characterize  par¬ 
ticular  features  of  this  equilibrium  surface.  , 

In  nonlinear  mechanics  the  principal  tools  for  such  a  computational  analysis 
are  the  so-called  incremental  methods.  These  procedures  were  ('•jveloped  more  or 
less  independently  in  the  engineering  literature.  But  they  are  now  also  re¬ 
cognized  to  be  closely  related  to  the  continuation  methods  used  for  some  time 
in  mathematics  in  general  and  in  numerical  analysis  in  particular.  The  liter¬ 
ature  in  this  area  is  extensive-  we  refer  only  to  [Zl]  for  a  discussion  about 
the  connection  between  incrcnk. .  approaches  for  structural  problems  and  con¬ 
tinuation  methods,  to  [ 8 ]  for  a  historical  overview  of  uses  of  continuation 
techniques  in  mathematics  and  to  [2],  [35]  for  some  literature  survey  of 
numerical  aspects  of  continuation  methods. 

Not  surprisingly  there  are  differences  between  the  methods  used  in 
structural  engineering  and  numerical  analysis  and  neither  is  directly  suited 
to  the  analysis  of  an  equilibrium  surface.  In  the  numerical  analysis  liter- 
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ature  continuation  methods  are  usually  considered  only  as  tools  for  de¬ 
termining  a  specific  solution  y*  of  a  given  nonlinear  operator  equation 
Gy  =  0.  For  this  the  equation  is  imbedded  into  a  one-parameter  family 
H(y,t)  =  0  which  has  a  solution  y  =  y(t)  for  each  fixed  t  in  some 
interval,  say,  0  _<  t  _<  1.  (See  eg.  [15]  for  a  survey  of  such  imbeddings.) 

If  y(t)  depends  continuously  on  t  and  satisfies  y(0)  =  y°  and  y(l)  =  y*, 
where  y°  is  a  known  point,  then  the  numerical  process  constructs  a  sequence 
of  points  in  the  proximity  of  the  path  y(t),  0^  t^  1,  starting  at  y° 
and  ending  at  the  desired  point  y*.  On  the  other  hand,  in  structural  me¬ 
chanics  incremental  methods  usually  are  designed  to  follow  numerically  a 
specific  load  curve  parametrized  by  a  load  intensity.  Hence,  while  in  the 
imbedding  approach  the  parameter  is  essentially  artificial,  in  the  incremental 
procedures  it  has  an  intrinsic  meaning  for  the  application,  and,  even  more 
importantly,  there  is  no  longer  a  fixed  endpoint  which  is  the  aim  of  the 
computation,  but  the  load  curve  itself  is  of  interest. 

For  a  numerical  analysis  of  a  given  equilibrium  surface  we  need  to  con¬ 
sider  continuation-methods  in  a  broader  sense  as  a  collection  of  numerical 
procedures  for  completing  at  least  the  following  three  basic  tasks: 

(i)  Follow  numerically  any  curve  on  the  surface  specified  by  a 
particular  combination  of  parameter  values  with  one  degree 
of  freedom. 

(1.1)  (ii)  On  any  such  curve  determine  the  exact  location  of  target 
points  where  a  given  state  variable  has  a  specified  value. 

(iii)  On  such  a  curve  identify  and  compute  exactly  the  critical 
points  where  stability  may  be  lost. 

Beyond  this  various  more  special  tasks  may  arise  as,  for  example,  the  following 


ones: 
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(iv)  From  any  one  of  the  critical  points  determined  under  (iii) 
follow  a  path  in  the  critical  boundary. 

(1  2) 

(v)  On  any  one  of  the  curves  (i)  determine  the  location  of 

bifurcation  points  and  the  paths  intersecting  at  that  point. 

Methods  which  are  either  directly  applicable  or  can  be  readily  adapted  to  com¬ 
pleting  these  various  tasks  have  been  proposed  by  various  authors.  In  particular, 
for  (i)  the  literature  is  very  large  and  we  refer  here  only  to  the  mentioned 
surveys  [2],  [35].  Methods  relating  to  (iii)  were  described,  for  instance,  in 
[1  1.  [20],  [22],  [33],  [34],  and  for  (iv)  and  (v)  we  refer  to  [28]  and  [10], 
[25],  respectively,  where  also  further  references  are  given. 

So  far  only  a  few  library  programs  for  performing  these  various  tasks  have 
been  published.  Without  claim  for  conpleteness  we  mention  here  [14],  [38]. 

Each  one  of  these  programs  has  the  objective  of  computing  a  specified  solution 
curve  of  a  nonlinear  equation  by  a  continuation  approach  along  the  lines  sketched 
above.  In  this  paper,  we  present  a  new  library  package  specifically  written  with 
the  objective  of  completing  the  three  basic  tasks  (1.1)  (i),  (ii),  (iii).  The 
package  can  be  expanded  to  incorporate  facilities  for  (1.2)  (iv),  (v),  but  this 
will  not  be  addressed  here.  The  package  is  based  on  the  continuation  approaches 
introduced  in  [26],  [27]  and  incorporates  some  of  the  concepts  of  steplength 
determination  discussed  in  [7].  At  the  same  time,  new  techniques  of  parameter 
adaptation  are  utilized  here  based  on  a  prediction  of  changes  in  the  curvature 
of  the  continuation  path. 

As  with  all  programming  packages  further  improvements  are  possible.  For 
example,  it  is  planned  to  introduce  an  automatic  first  step  selection  and  a 
function-scaling  option.  Special  versions  incorporating  facilities  for  the  tasks 
(1.2)  are  also  being  designed.  But  since  all  these  changes  are  built  on  the  pre¬ 
sent  package  the  presentation  of  a  documentation  of  PITCON  in  its  basic  form 
appeared  desirable  and  justified. 


\ 
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BASIC  FORMULATION 


Generally,  after  suitable  discretizations,  the  equilibrium  problems 
mentioned  in  the  introduction  lead  to  a  finite-dimensional,  non-linear 
equation  of  the  form 


(2.1)  G{y,p)  =  0 

where  y  e  R  is  a  vector  of  state  variables,  p  e  R  a  vector  of  parameters, 
and  G:  R*"  x  -*■  R*"  a  given  function.  Then  we  are  interested  in  the  features 
of  the  set 


(2.2)  E(G)  =  {(y,p)  e  r"  x  R»*;  G(y,p)  =  0} 

of  all  solutions  of  (2.1).  Under  well-known  conditions  E(G)  represents  an 
r-dimensional  manifold  in  r"*  x  r'*. 

In  most  applications,  interest  centers  on  tracing  paths  on  E(G)  which 
are  characterized  by  r-1  relations  between  the  parameters.  In  other  words, 
we  are  given  a  suitable  mapping  K:  R*^  -*■  R*^"^  and  wish  to  compute  the  subset 
of  E(G)  defined  by  the  augmented  equations 

G(y,p)  =  0, 

(2.3) 

Kp  =  0. 

In  this  formulation  we  should  include  the  parameters  in  the  list  of  variables, 
in  which  case, (2.3)  represents  a  system  with  one  more  variable  then  equations. 
Then,  for  ease  of  notation,  it  is  reasonable  to  combine  the  vectors  y  and 
p  into  one  vector  x  of  dimension  n  »  m+r.  Moreover,  from  the  viewpoint 
of  our  package  of  library  programs  it  is  natural  to  assume  that  both  mappings 
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G  and  K  of  (2.3)  are  provided  for  by  the  user.  In  other  words,  we  may 
write  (2.3)  as  one  equation 

(2.4)  Fx  *  0 

with  a  user-specified  mapping  F;  r"  r”~^.  Note,  however,  that  in  this 
underdetermined  equation  (2.4)  no  one  variable  is  explicitly  identified  as 
continuation  variable  as  is  typical  in  the  incremental  and  continuation 
methods  mentioned  in  the  introduction. 

We  assume  here  that  the  given  mapping  F  has  the  following  properties: 

(i)  F  is  continuously  differentiable  on  R*’. 

(2.5)  (ii)  The  derivative  DF(x)  of  F  is  locally  lipschitzian  on  r". 

(iii)  The  regularity  set  R(F)  =  {x  e  r"^;  rank  DF(x)  *  n-1}  is 
non-empty  and  therefore  an  open  subset  of  r”. 

From  (2.5)  it  follows  (see  [26])  that  the  tangent  map  specified  by 

n  /DF(x)\ 

(2.6)  T:  R(F)  ^  r",  DF(x)Tx  =  0,  l|TxiL  =  l,  det  (  ,  >0 

^  UTx)V 

is  uniquely  determined  and  locally  lipschitzian  on  R(F).  Furthermore,  (2.5) 
implies  that  the  regular  solution  set  6(F)  O  R(F)  of  F  is  either  empty 
or  a  one-dimensional  C^-manifold  in  the  open  set  R(F).  Our  objective  is  to 
determine  numerically  a  non-empty  connected  component  E*  of  E(F)  n  R(F). 

It  is  well-known  (see  eg.  [17])  that  such  a  component  E*  is  diffeomorphic 
either  to  the  circle  or  to  some  interval  (that  is,  some  connected  subset)  of 
r\  Hence,  E*  is  uniquely  determined  by  any  one  of  its  points  x°  e  E(F)  D  R(F) 
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and  we  denote  this  by  writing  E*(F,x°)..  Note  that  for  any  e  E*(F,x°) 
we  have  E*(F,x^)  =  E*(F,x°). 

A  parametrization  by  arclength  of  E*(F,x°)  is  a  solution  of  the 
initial  value  problem 

(2.7)  X  =  Tx,  x^O)  =  x°. 

Note  that,  since  T  is  locally  Lipschitzian,  (2.7)  has  a  unique  solution 
which  cannot  terminate  inside  R(F).  Evidently  standard  ODE-solvers  may 
be  applied  to  solve  (2.7)  numerically.  This  has  been  pursued  for  some 
time  in  the  literature  (see  eg.  [3],  [el.  [131.  [37]).  Independent  of 
this, the  choice  of  the  arclength  for  the  parametrization  of  E*(F,x°)  has 
been  proposed  by  many  authors.  Notably  H.  B.  Keller  and  his  co-workers  (see 
eg.  [lOl,  [111)  have  advocated  this  choice  for  some  time.  It  is  also 
the  basis  of  incremental  procedures  given  in  [5l,  [29]  and  has  been  more  or 
less  implicit  in  various  papers  in  the  field. 

Our  programs  here  are  based  more  generally  on  the  structure  of  E*(F,x°) 
as  a  one-dimensional  manifold  and  use  a  local  parametrization  at  each  point 
computed  along  E*(F,x°).  A  natural  class  of  such  local  parameters  are  the 
n  components  of  the  vector  x.  We  call  a  process  based  on  this  choice  of 
parametrization  a  locally- parametrized  continuation  method. 
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3.  OUTLINE  OF  THE  PROCESS  AND  BASIC  STEPS 

As  noted  before,  our  objective  is  to  determine  numerically  a  non-empty 
component  E*{F,x°)  of  the  regular  solution  set  B{F)  R(F).  For  the 
discussion  it  is  useful  to  consider  a  parametrization  by  arclength  of 
E*(F,x°),  that  is,  a  function  x:  J  ->■  E*(F,x°)  which  maps  some  interval 
jc  diffeomorphically  onto  some  open  subset  of  E*(F,x°)  such  that 
llx(s)ll2  *  1  for  s  e  J.  We  may  assume  also  that  x(0)  =  x°,  0  e  J. 

The  process  described  here  belongs  to  the  class  of  predictor- corrector 
continuation  methods.  Starting  from  x°  it  produces  a  sequence  of  approxi- 
nations  x  «  x(S|^),  k  =  0,1,...,  corresponding  to  some  sequence 
0  =  $0  <  Si  <  Sg  <  . . .  of  arclength  values.  Note,  however,  that  in  general 
the  values  s^,S2,..-  are  only  approximately  computable  and  are  of  limited 
interest  in  most  applications. 

In  our  program  the  principal  steps  performed  during  one  continuation 
step  are  as  follows: 

1.  Initialization. 

2.  Check  for  and  computation  of  target  point,  if  desired. 

3.  Calculation  of  tangent  vector  and  determination  of  new  local 
continuation  parameter. 

(3.1)  4.  Check  for  and  computation  of  limit  point,  if  desired. 

5.  Steplength  computation. 

6.  Computation  of  predicted  point  and  corrector  iteration. 

7.  Storage  of  data  and  return. 


The  sequencing  of  these  steps  is  dictated  by  the  data-flow.  For  the  de¬ 
scription  of  the  details  it  will  be  advantageous  not  to  adhere  to  this 
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sequence.  Instead,  in  the  remainder  of  this  section,  we  discuss  the  basic 
steps  3.  and  6.  Then  the  next  section  introduces  the  new  steplength  compu¬ 
tation  used  in  step  5.  and  section  5  covers  steps  2.  and  4.  The  data- 
handling  steps  1.  and  7.  should  be  self-explanatory  from  the  documentation 
of  the  program  itself. 

Let  e\...,e'^  be  the  natural  basis  vectors  of  Then  it  is  readily 
verified  that  (see  eg.  [26]) 


(3.2) 


=  [(e^)^Tx]  det 


V  X  e  R(F) ,  i  =  1 ,. . . ,n. 


where  the  matrix  occurring  on  the  right  is  non-singular.  Hence,  for  any  in¬ 
dex  i,  1  1  i  ^  n,  such  that  (e’)^Tx  /  0,  the  solution  v  e  r'^  of  the 
linear  system 


(3.3) 

is  uniquely  defined. 


Evidently,  then 


(3.4) 


» 


and,  in  line  with  (2.6),  we  should  set 
(3.5)  a  =  sign(v^e^)  sign  det 

i 


As  long  as  the  solution  path  remains  completely  in  R(F)  this  is  satisfactory. 
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But  frequently  in  applications  we  may  encounter  a  bifurcation  point  x*  ^  R(F) 
where  several  solution  paths  terminate.  For  example,  using  arclength  repre¬ 
sentations  we  may  find  solutions  x"^:  J.c  -*-R(F),  j  =  1,...,4,  for 
which  x'^(s)  tends  to  x*  when  s  tends  to  one  of  the  endpoints  of  J.. 

Moreover,  it  often  happens  that  there  are  pairs  of  these  solutions,  say,  x^ 

2  12 
and  X  for  which  lim  x  (s)  =  -lim  x  (s) 

at  X*,  (see  Fig.  1).  In  other  words,  if 

we  disregard  the  direction  of  the  solutions, 

they  appear  to  form  one  smooth  curve  through 

X*.  In  such  a  case,  when  the  process  moves 

along  x^  toward  x*  it  usually  "jumps" 

2 

over  X*  onto  x  .  Then,  unless  we  reverse 
the  sign  of  a  in  (3.4)  the  tangent  will 
again  point  toward  x*  and  the  process  reverses  direction. 

In  order  to  avoid  this  problem  suppose  that  the  point  x  in  (3.4)  is  the 
k-th  approximation  computed  along  the  curve.  Then  a  is  determined  as  follows 

I  dir,  if  k  =  0 

+1  ,  if  sign  v^e^  =  sign(Tx'^"^)^e’ 

-1  ,  otherwise 

where  dir  is  a  user  specified  direction  at  the  starting  point.  By  comparing 
this  value  of  a  with  that  of  (3.5)  we  can  detect  if  the  process  did  jump 
over  a  bifurcation  point  of  odd  multiplicity.  Obviously,  bifurcation  points 
of  even  multiplicity  cannot  be  found  this  way. 

Once  the  tangent  Tx  has  been  obtained  we  determine  the  indices  j-j 
and  ^2  largest  and  second  largest  component  of  Tx  in  modulus. 


Figure  1 
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respectively.  The  relation  (3.2)  certainly  suggests  that  the  index  ij^, 

^  f>,  of  the  new  local  continuation  variable  be  set  equal  to  j,| . 

However,  if  we  are  approaching  a  limit  point  in  the  j-j-th  variable  then 

this  choice  may  be  disadvantageous.  Accordingly,  if  the  following  three 

conditions  are  simultaneously  satisfied 


(i) 


$ 


(3.7)  (ii)  !(e^2)T^^kj  ^  ^ 


(iii)  |(e  2)'^Tx‘^|  >  y|(e'')' Tx'^l. 


1^T.  k, 


with  a  fixed  y,  0  <  y  <  1,  then  we  set  ij^  =  Of  course,  if  we  don't 

have  a  previous  tangent  vector  this  check  has  to  be  bypassed.  The  new  con- 

k+1 

tinuation  index  ij^  will  be  used  for  the  computation  of  the  next  point  x 

k+1  n 

and  its  tangent  Tx  .  For  the  tangent  computation  at  x  a  continuation 

index  is  assumed  to  be  given  by  the  user. 

With  the  tangent  Tx  and  the  steplength  h|^  >  0  determined  by  the 

steplength  algorithm  of  section  4  we  compute  now  the  predicted  point 
^  k  k  k 

X  =  X  +  h|^Tx  .  Then  any  appropriate  iterative  method  for  the  solution  of 
the  augmented  equation 


(3.8) 


Fx  = 


Fx 


(e''')^(x-x''), 


=  0 


^  k 

starting  from  x  may  be  used  as  a  corrector  process.  In  the  program  we 
use  either  the  regular  Newton  method  or  its  modified  form  in  which  the 
Jacobian  at  the  starting  point  is  held  fixed. 


n 


0  k  2 

Let  y  =x,y‘,y,...  be  the  iterates  produced  in  this  way.  The 
process  has  to  incorporate  provisions  for  monitoring  the  convergence  and 
for  aborting  the  iteration  as  soon  as  divergence  is  suspected.  In  the 
program  non-convergence  is  declared  if  any  one  of  the  following  three 
conditions  is  true 

(i)  1 1  Fy'^l  I  i  6|  I  I  for  some  j^l, 

(3.9)  (ii)  I  I  i  e|  I  tor  some  j  >  2, 

(tH) 

For  the  constant  6  we  use  9  =1.05  except  in  the  first  check  of  (3.9)  (i) 

where  9=2  is  chosen.  The  maximal  iteration  count  depends  on  the 

max 

method.  For  the  regular  Newton  process  we  set  =  10  and  double  this 

for  the  modified  method.  In  the  case  of  non-convergence  the  predictor  step 
is  reduced  by  a  given  factor,  for  example  1/3,  unless  the  resulting  step  is 
below  a  given  minimal  steplength. 

Convergence  is  declared  if  either  one  of  the  two  conditions  holds  for 
an  iterate: 

(i)  I  |FyJ|  I  <  8  for  some  j  >  0, 

(3.10) 

(ii)  (I  iFy^M  1  e^jjj)  and  (|  ly^V^ 1 1  1  E^bs  ^rel ' I  ^ 

The  tolerances  Cgbs’  ^rel  specified  and  is  the  smallest 

floating  point  number  such  that  1.  =  1.  +  ^niach’  tests  (3.9)  and 

(3.10)  the  maximum  norm  is  used. 


4.  THE  STEPLENGTH  ALGORITHM 


For  the  points  x  ,  k  =  0,1,...,  approximating  the  continuation  curve 
x:  J  -►  E*(F,x  )  the  achievable  error  |lx  -x(S|^)||  is  solely  determined  by 
the  termination  criterion  (3.10)  of  the  corrector  process.  In  contrast  to 
this  the  standard  ODE-solvers  involve  a  corrector  equation  obtained  by  extra¬ 
polation  for  which  the  solutions  are  not,  in  general,  on  the  exact  curve.  As 
a  consequence  the  available  error  for  the  ODE-solvers  depends  on  the  history 
of  the  process  up  to  that  point,  and  this  in  turn  has  a  strong  influence  on 
the  step- select! on.  On  the  other  hand,  for  our  continuation  process  any  step 
h|^  >  0  along  the  Euler  line  is  acceptable  in  principle  if  only  the  corrector 
converges  from  the  predicted  point  x*'.  Moreover,  in  [261.it  was  shown  that 
any  compact  segment  of  the  continuation  curve  in  R(F)  has  an  e-neighborhood 
for  some  e  >  0  in  which  Newton's  method  will  converge  to  the  curve. 

This  suggests  that  we  estimate  the  radius  of  convergence  of  the  corrector 
process  at  the  computed  points  and  extrapolate  these  radii  to  the  next  point 
about  to  be  determined.  In  practice  the  estimate  of  a  convergence  radius  at 
some  continuation  point  would  have  to  be  based  on  the  corrector  iterates  which 
led  to  that  point.  Unfortunately,  as  was  proved  in  [7],  this  represents  in¬ 
sufficient  information  for  obtaining  such  an  estimate.  On  the  other  hand,  an 
approach  was  presented  in  [7  ]  which  allows  for  an  assessment  of  the  con¬ 
vergence  quality  of  the  particular  sequence  of  corrector  iterates. 

For  details  of  this  approach  we  refer  to  the  cited  article.  In  brief, 
let  {y^}  be  a  given  sequence  with  limit  y*  generated  by  an  iterative 
process  and  denote  the  errors  by  e^  =  !|y^-y*lli  i  =0,1,...  .  The 
definition  of  any  convergence  measure  is  based  on  a  hypothetical  model  of 
the  behavior  of  the  errors.  For  example,  if  (y’}  converges  linearly  it  is 
reasonable  to  assume  that 
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(4.1)  0  j<  <  Xe^. ,  i=0,l,... 

with  some  constant  x,  0  <  X  <  1,  depending  on  {y^}.  Suppose  now  that 

i  * 

the  process  was  terminated  with  the  iterate  y  .  Then 


(4.2) 


i* 


>  2, 


represents  a  computable  estimate  of  x. 

In  the  setting  of  our  continuation  process  suppose  now  that  the  y\ 
i  =0,1,...,  are  the  corrector  iterates  leading  from  the  current  predicted 

Ak  0  k+1  1^* 

point  X  =  y  to  the  new  continuation  point  x  =  y  .  Then 


(4.3) 


4^  - 


is  the  correction-distance.  For  the  modified  Newton  method  the  convergence 
is  indeed  linear,  and  a  reasonable  aim  in  the  construction  of  the  steps  along 
the  curve  is  to  ensure  that  the  number  of  corrector  iterates  remains  about 
constant.  In  other  words,  we  aim  at  taking  always,  say,  m*  corrector  steps. 
Hence,  under  the  heuristic  assumption  that  the  error  model  (4.1)  remains 
valid  for  some  interval  of  starting  errors  e^^  around  6|^,  we  should  have 
begun  with  an  "ideal  starting  error"  6j^  =  such  that 


(4.4) 


-X"*  4.  4  ,1*  4^ 


and  therefore 
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(4.5) 


-  A  =  OJ 


In  our  program  we  use  m*  =  10  for  the  modified  Newton  method  and  enforce 
always  that  0.125^  1  8. 

This  technique  is  also  readily  applicable  for  Newton's  method.  In  [ 7  ] 
two  different  hypothetical  error  models  for  the  Newton  process  were  discussed. 
Here  we  use  only  one  of  these  models,  namely  the  one  arising  in  the  attraction 
theorem  formulated  in  [24].  In  essence,  under  certain  conditions 
about  the  equation  and  the  desired  limit  y*  of  the  Newton  process  there 
exists  a  radius  r*  >  0  such  that  for  any  starting  point  y°  in  the  ball 
B(y*,r*)  the  relative  errors  e.  =  e^./r*,  i  =  0,1,...,  satisfy 


(4.6)  0  ^  S-j+l  ^  'l>(s.j ) »  i  =  0,1 .... ,  (J)(t)  =  ■jrjt  ♦  0  £  t  ^  1 . 


The  radius  r*  depends  on  global  information  about  the  equation  and  is  not 
accessible.  If  0  <  <  1  and  the  {e^.}  satisfy  (4.6)  then  we  have 


(4.7)  e.  < 

’  “  ^  °  1+2  cosh  2a 


where  a  is  the  unique  positive  solution  of 


,  i  =  0,1,...,  Hq  =  Ep, 


(4.8) 


tp(a)  =  Hq.  <i'(ct)  = 


Moreover,  for  any  m,  0  <  o)  <  1 ,  and  i*  ^  2  the  equation 
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(4.9) 


1 

'i'(a) 


1  +  2  cosh  a 
- ^ 

1+2  cosh  2^  a 


has  a  unique  solution  a  >  0. 

Now  suppose  that  {yM  denotes  the  sequence  of  Newton  iterates  and 
that  the  process  was  terminated  at  y  .  As  in  the  linear  case  we  use  the 
approximati on 


(4.10) 


lly'  -y°ll 


^i*-l 


-  n. 


and  compute  with  this  oi  the  solution  a  of  (4.9)  which  gives  the  estimate 
Hq  =  'j^(a)  of  Eq.  Now  we  proceed  as  before  and  obtain  the  factor 


(4.11) 


for  the  ideal  starting  error  by  determining  the  unique  solution  fi^,  0  <  <  1, 

of 


(4.12)  ^*(0;)  =  <D’*(np). 

Since  the  iterates  are  explicitly  known  the  various  equations  are 
not  difficult  to  solve  numerically.  However,  for  the  computation  it  is  more 
advantageous  to  introduce  a  least  squares  fit  of  as  a  function  of  us 
for  all  relevant  values  of  i*.  In  the  program  we  use  m*  =  4  and  the 
approximations  for  9(^  given  in  Table  1.  Note  that  as  before  we  restrict 
to  the  interval  0.125  ^  £  8. 


0)  e 


a 

b 

K 

0.8735115 

i 

1 

i 

1 

0.1531947 

i  0.8735115 

0.9043128  -  0.7075675  In  oo 

2 

0.03191815 

1 

j  0.1531947 

-4.667383  -  3.677482  In  u 

0 

0.03191815 

8 

0.4677788 

1 

1 

0.6970123(-3) 

0.4677788 

0.8516099  -  0.1953119  In  oi 

3 

1 

0.1980863(-5) 

0.6970123(-3) 

-4.830636  -  0.9770528  In  oj 

0 

0.1980863(-5) 

8 

4 

0 

1 

1 

0.3339946(-10) 

1 

1.040061  +  0.03793395  In  uj 

5 

0 

0.3339946(-10) 

0.125 

0.1122789(-8) 

1 

1.042177  +  0.04450706  In  0) 

6 

0 

0.1122789(-8) 

0.125 

>  7 

0 

1 

0.125 

Table  1 


We  turn  now  to  the  algorithm  for  the  determination  of  the  steplength 

Ic  k 

h|^  >  0  along  the  Euler  line  ir(t)  =  x  +  t  Tx  used  for  the  prediction.  In 
order  to  estimate  the  distance  between  7r(t)  and  the  exact  curve  x  =  x(s) 
we  introduce  the  quadratic  Hermite-Birkhoff  interpolation  polynomial 
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(4.13) 


q(t)  =  +  t  Tx'^  +  ^  w'J.  ^  (Tx'^-Tx'^’b,  As,^  =  |  |  1 2 


for  which 


(4.14) 


q(0)  =  x*".  q'(0)  =  Tx^  q'(-AS^)  =  Tx''"^ 


Since 


(4.15) 


k 


w 


rl 

x"(S|^-aAS|^)do  =  x"(S|^-oAS|^), 

0 


0  <  o  <  1 , 


the  quantity 


(4.16)  Il«'^ll2  *  a!;;’  I“kl’  Ok  =  arccos  ({Tx'^)'^Tx'^'‘') 

represents  an  approximation  of  the  curvature  of  the  exact  point  at  some  point 
between  x(S|^  .|)  and  x(S|^). 

It  is  tempting  to  derive  from  q  a  prediction  of  the  curvature  to  be 

expected  during  the  next  continuation  step.  However,  a  closer  computation 

k  2  1 

shows  that  the  value  of  the  curvature  of  q  assumes  its  maximum  !|w  II2/COS  ^ 

at  t  AS|^  and  that  for  increasing  t  this  value  decreases  rapidly.  For 

k  1 

example,  at  t  =  0  the  curvature  of  q  equals  only  |  |w  Ijg  |cos  j  and 

for  positive  t  no  reasonable  predictive  information  can  be  gained  this  way. 

The  relation  (4.15)  suggests  the  use  of  the  simple  linear  extrapolation 


tent  _ 
^k  “ 


AS. 


‘h  *  “k-i 


(I  I""] 


2 


k-1 


I.) 


(4.17a) 


for  a  prediction  of  the  curvature  during  the  next  continuation  step.  However, 
this  value  may  become  negative  and  accordingly  we  use  instead 


(4.17b)  y,  .  max  |y„,„.y*"*) 

with  a  given  small  ^  0* 

Most  of  the  data  discussed  so  far  are  sketched  in  Figure  2.  In  order 
to  derive  a  formula  for  the  desired  predictor  step  hj,  we  note  that 

(4.18)  llq(t)  -  TT(t)||2  =  ^  t2l|w''||2 

represents  an  estimate  of  the  distance  between  the  Euler  line  and  the  exact 
curve.  In  fact,  for  smooth  curves  the  error  of  this  estimate  is  asymptotically 
or  order  three  in  max  ( |  tj  ,aS|^)  as  this  quantity  tends  to  zero.  Hence,  if  we 


want  this  distance  to  be  at  most  equal  to  a  tolerance  e  >  0  then  we  should 

choose  the  next  step  as 

(4.19) 

— 

/ 1  Iw  1 

'2 

It  is  natural  to  replace  the  curvature 

|jw  ((2  by  the  predicted  value  Y|^ 

of  (4.17)  and  to  relate  the 

tol erance  e 

to  the  "ideal  starting  error"  6j^ 

obtained  earlier.  As  Figure  2  indicates 

it  is  unreasonable  to  expect 

Hence,  we  use  instead 

» 

^min  ^^k 

(4.20)  =  ( 

if  8j;>bs„ 

sa 

otherwise 
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with  a  small  >  0,  e.g.,  =  0.01.  Then  a  tentative  predicted  step 

is  given  by 


(4.21) 


From  the  form  (3.8)  of  the  augmented  equation  we  see  that  the  corrector 
iterates  remain  in  a  hyperplane  perpendicular  to  the  basis  vector  e  through 
the  predicted  point.  Then  Figure  2  suggests  that  we  adjust  the  predicted 
steplength  hj^  so  as  to  ensure  that  hj^  will  be  approximately  equal  to 
There  is  no  need  to  enforce  this  too  rigidly.  It  suffices  to  define  a  new 
tentative  step  by  the  requirement 


(e’V  ir(h^2))  ,  (gijT  ^(^(1)^ 


whence. 


(4.22) 


2as. 


(1 


Tx 


k-1 


(e’)^ 


Tx' 


-)1. 


This  formula  may  involve  subtractice  cancellation  and  has  to  be  evaluated  in 
double  precision. 

f  21 

The  final  value  hj^  of  the  steplength  is  now  obtained  from  hj^  '  by 

enforcing  three  different  bounding  requirements.  First  of  all,  if  the  pre- 

k-1  k 

vious  continuation  step  from  x  to  x  was  obtained  only  after  a  failure 
of  the  correction  process  and  a  corresponding  reduction  of  the  predicted  step, 
then  we  should  not  allow  hj^  to  exceed  aS|^.  Secondly,  as  in  the  ODE  solvers 
we  need  to  control  both  the  relative  growth  and  the  absolute  size  of  the  pre- 


dictor  step.  Thus,  we  require  that 


(4.23) 


h^  <  <AS^, 


h  .  <  h,  <  h 

mm  —  k  —  max 


where  k  is  some  factor,  say,  <  =  3,  and  h  h_^^  depend  on  the  machine 

min  max 

as  well  as  the  requirements  of  the  problem.  It  should  be  obvious  how  the 

(2) 

final  step  hj^  is  obtained  from  hj^  '  on  the  basis  of  these  restrictions. 


22 


5.  THE  COMPUTATION  OF  TARGET  AND  LIMIT  POINTS 

By  generating  a  sequence  of  solution  points  on  a  given  curve,  the  con¬ 
tinuation  process  reveals  the  shape  of  the  curve,  but  there  are  often  other 
items  of  interest  that  need  to  be  studied  as  well.  Our  program  is  designed 
to  pause  during  the  continuation  steps  in  order  to  seek  out  special  points 
that  the  user  has  requested,  namely,  target  and  limit  points. 

A  target  point  x  £  E*(F,x*^)  is  a  point  on  the  solution  curve  for 
which  the  component  x^.  =  (e^)^x  with  given  index  i  =  IT  has  a  prescribed 
value  x^.  =  XIT.  Limit  or  turning  points  with  respect  to  a  given  index 
i  =  LIM  are  points  x  e  E*(F,x®)  where  the  i-th  component  (e’)^Tx  is  zero. 
More  specifically,  since  it  is  computationally  unreasonable  to  attempt  to 

compute  zeroes  of  even  order,  we  are  concerned  only  with  limit  points  on  the 

i  T 

continuation  curve  where  (e  )  Tx{s)  changes  sign. 

It  might  be  mentioned  that  bifurcation  points  represent  another  inter¬ 
esting,  special  class  of  points.  But  in  that  case  we  are  not  only  interested 
in  the  specific  location  of  the  point  but  also  in  the  solution  curves  tha\, 
branch  off  from  it.  This  is  exactly  the  task  (1.2)  (v)  listed  earlier.  The 
corresponding  procedures  (loc.  cit.)  would  add  considerably  to  the  complexity 
of  our  program,  and,  since  their  utility  tends  to  be  of  a  more  specialized 
nature,  it  was  decided  not  to  cover  task  (1.2)  (v)  (nor  (1.2)  (iv))  in  the 
present  program. 

As  indicated  before,  the  determination  of  a  target  or  limit  point  re¬ 
presents  an  interruption  in  the  normal  flow  of  the  continuation  program. 

After  at  least  one  step  has  been  taken,  the  program  has  available  an  old 

k-1  k  k-1 

point  X  ,  a  new  point  x  and  the  tangent  vector  Tx  .  Normally,  then 

we  turn  to  the  computation  of  Tx  ,  of  the  new  steplength,  and  finally  of  the 
k+1 

next  point  x  .  But  if  the  index  IT  or  LIM  is  non-zero  then  these 


computations  are  postponed  for  the  search  of  a  target  or  limit  point, 
respectively.  We  discuss  these  cases  separately: 

Target  points:  Suppose  that  a  non-zero  value  of  i  =  IT  and  associated 

i  T  k-1 

value  =  XIT  have  been  given.  If  lies  between  (e  )  x  and 

1  T  k  0 

(e  )  X  ,  then  it  is  assumed  that  a  solution  point  x  e  E*(F,x°l  with 

(e^)^x  =  x^.  is  nearby.  In  this  case  a  point 

(5.1)  y{t)  =  {l-t)x''‘’  +  tx'^.  0<t<l, 

k-1  k  i  T  - 

on  the  secant  between  x  and  x  is  determined  such  that  (e  )  y(t)  =  x. . 

Now  with  the  augmenting  equation  (e’)^x  =  x^.  the  corrector  process  is  applied, 

and,  if  it  terminates  successfully  the  resulting  point  is  taken  as  the  desired 

target.  Otherwise,  a  failure  is  indicated  for  the  target  routine.  In  either 

case,  the  routine  returns  and  on  the  next  call  the  continuation  loop  will  pick 

up  from  where  it  was  interrupted.  Note  that  in  effect  the  target  routine  uses 

k-1 

the  IT-th  variable  for  the  local  parametrization  of  the  curve  between  x 
and  X  .  This  may  be  an  inferior  choice  of  parameter  for  the  corrector  but 
it  allows  us  to  enforce  that  the  resulting  target  point  x  e  E*(F,x°)  indeed 
satisfies  (e^)^x  =  x. .  Clearly,  for  very  large  continuation  steps  we  have 
no  guarantee  that  all  target  points  will  be  detected  or  that  a  target  compu¬ 
tation  will  succeed.  Thus,  the  utility  of  the  target  routine  will  depend  on 
the  maximal  allowed  stepsize  that  has  been  chosen. 

Limit  points:  If  the  limit  point  index  i  =  LIM  is  non-zero,  then  a  limit 
point  determination  is  carried  out  after  a  target  point  search  has  been 
successfully  or  unsuccessfully  completed,  provided  it  was  called  for  at  all. 
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k-1  k 

Recall  that  we  still  have  as  current  information  the  vectors  x  ,  x  , 

k-1  k 

and  Tx  .  Now  the  new  tangent  Tx  is  evaluated  and  if 

sign(e’ i  sign(e’ Tx*^  for  i  =  LIM  0)  then  a  limit  point 

search  is  begun.  For  this  the  index  i  of  the  largest  component  in 

k  k-1 

modulus  of  the  secant  direction  x  -  x  is  chosen  as  a  local  para- 

k-1  k 

metrization  of  the  curve  between  x  and  x  .  More  specifically,  suppose 

that  x:  [S|^  represents  the  segment  of  the  curve  between 

k-1  k 

X  and  X  .  Then  i  is  assumed  to  be  the  index  of  a  local  coordinate 

for  which  there  exists  a  bijective  parameter  transformation  4):  [0,1]  -*■  [S|^_i,S|^] 

such  that  (e’Vy(t)  =  (e^)^  x(<))(t)) ,  0  _<  t  _<  1 ,  where  y(t)  is  defined  by 

(5.1) . 

Hence,  we  may  consider  the  function 

(5.2)  g:C0,ll^R\  g(t)  =  (e’ )’’’ Tx((^(t) ) ,  0<t<l. 

and  our  problem  is  to  determine  a  zero  of  g.  Since  by  assumption 

sign  g(0)  f  sign  g(l),  a  rootfinder  of  the  Dekker-Brent  type  can  be  applied. 

For  the  evaluation  of  g(t)  we  use  the  augmenting  equation  (e^)^x  =  (e’)^y(t) 
and  apply  the  corrector  process  with  y(t)  as  starting  point.  If  it  terminates 
successfully  with  some  x  then  Tx  can  be  evaluated  and  we  set  g(t)  =  (e’)^Tx. 
Hence,  g  is  certainly  costly  to  compute  and  we  require  an  efficient  root- 
finder  to  speed  the  conve'-gence  of  the  limit  point  routine.  A  specially 
modified  version  of  the  routine  given  in  [4]  is  used  in  our  program.  Clearly, 
as  in  the  case  of  target  points,  we  may  fail  to  detect  a  limit  point  if  the 
continuation  steps  are  too  large  and  in  such  a  situation  the  rootfinder  may 
also  fail  to  converge.  In  addition,  the  evaluation  of  g  may  run  into 
difficulties  when  the  desired  limit  point  is  near  a  bifurcation  point. 
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6.  SOME  NUMERICAL  EXAMPLES 

The  programs  described  here  have  been  used  extensively  with  excellent 
success  on  problems  from  many  different  areas.  We  include  here  only  a  few 
numerical  examples  to  illustrate  the  operation  of  the  programs. 


Exampl e  1 .  In  order  to  present  some  details  of  the  performance  of  the  pro¬ 
grams  we  consider  first  a  very  small  problem  which  was  originally  formulated 
in  [9  ]  and  subsequently  used  as  a  test  case  by  many  authors.  The  mapping  F 
has  here  the  form 


(6.1) 


Fx  = 


/ 


Xi  -  ^2  +  5x2 


2x2  +  34x2  -  47 


V  X  e  R' 


x.|  +  x|  +  x|  -  14x2  +  10x2  ' 


For  the  starting  point  x°  =  (15,-2, 0)^  the  s>'/^l!^l•')n  cv'-'^  passes  through 
X*  =  (5,4,1)^  and  this  point  is  chosen  as  target. 

Tables  2  and  3  show  runs  with  the  full  Newton  method  and  modified  Newton 
method,  respectively  as  corrector  process.  A  starting  step  h^  =  0.3  and 
maximum  step  =  25.0  were  used.  The  performance  for  the  two  correctors 

is  practically  the  same  although  the  step-prediction  exhibits  certain  differences 
due  to  our  assessment  of  the  corrector  distance.  Clearly,  the  use  of  the 
modified  Newton  process  is  much  less  expensive  and  hence  preferable  as  Table 
4  shows  which  sumnarizes  the  total  number  of  function  and  Jacobian  calls  in¬ 
cluding  those  for  the  target  calculation.  Comparative  performance  data  given 
in  [ 7 ]  for  this  problem  involved  22  continuation  steps,  15  step  reductions 
and  128  Jacobian  evaluations.  The  procedure  discussed  in  [191  required  25 
continuation  steps  but  no  further  details  were  provided  in  the  paper. 
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Continuation  point 

Total 

Conti n. 

Correct. 

^2 

^^3 

Variable 

Steps 

Comments 

15.000 

-2.00000 

0.00000 

^3 

— 

14.705 

-1 .9421 

0.065381 

’'l 

2 

14.285 

-1.7291 

0.26874 

^3 

3 

16.906 

-1.2094 

0.54684 

^2 

2 

24.918 

-0.59908 

0.55514 

’‘l 

3 

1 

48.974 

0.71803 

-0.080758 

’‘l 

3 

57.928 

1.2846 

-0.40736 

4 

Step  red. 

60.052 

1.5709 

-0.54035 

^1 

4 

Step  red. 

61.666 

2.0010 

-0.66683 

"2 

2 

-5.1039 

4.1510 

1.3464 

X2 

2 

Target  passed 

TOTAL 

25 

Computation  of  target 

4 

Table  2 


15.000 

-2.0000 

0.00000 

"3 

14.710 

-1 .9421 

0.065381 

"l 

14.285 

-1 .7291 

0.26874 

^3 

16.906 

-1.2094 

0.54685 

X2 

24.918  i 

-0.59906 

0.55514 

’^l 

48.975 

0.71810 

-0.080804 

"1 

57.289 

1.2847 

-0.40742 

"1 

60.053 

1.5711 

-0.54042 

!  ^1 

1 

61.666 

2.0013 

-0.66689 

Xg 

-4.4239 

4.1413 

1 .3229 

X2 

TOTAL 

Computation  of  target 
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Corrector  Process 

Newton 

Mod.  Newton 

Function  calls 

41 

53 

Jacobian  calls 

38 

21 

Table  4 


It  may  be  noted  that  the  solution  curve  has  two  limit  points  each  with 
respect  to  and  x^.  The  two  step  reductions  are  almost  unavoidable  here 
since  the  curve  has  a  long  straight  segment  followed  by  a  very  sharp  bend. 

The  target  computation  is  relatively  expensive  since  the  last  step  is  extremely 
large  due  to  another  straight  curve  segment. 

Example  2.  Maneuvering  airplanes,  especially  at  high  angles  of  attack,  some¬ 
times  undergo  sudden  jumps  in  their  response  to  the  pilot’s  control  inputs. 

The  problem  has  been  discussed  extensively  in  the  literature,  see,  for  example, 
[18],  [30],  [39].  Without  going  into  further  details  we  use  here  a  simplified 
version  of  a  system  of  five  equilibrium  equations  involving  the  roll  rate  (x^), 
pitch  rate  yaw  rate  (x^),  (incremental)  angle  of  attach  (x^),  side  slip 
angle  (xg),  elevator  angle  (xg),  aileron  angle  (xy),  and  rudder  angle  (xg). 
More  specifically,  for  the  particular  aircraft  discussed  in  [18]  these  equations 
have  the  dimensionless  form 

(6.1)  Fx  =  Ax  +  <j)(x)  =0,  V  X  e  R® 


where 


<l>(x)  = 


-0.727  X2X2  +  8.39  -  684.4  x^Xg  +  63.5  x^x^ 


0.949  x-jX^  +  0.173  x^Xg 


-0.716  x^X2  -  1.578  x^x^  +  1.132  x^x^ 


'  "l"2 


Figure  3  shows  some  solution  curves  on  the  three-dimensional  equilibrium  sur- 

O 

face  in  R  .  More  specifically,  in  all  cases  we  fixed  a  value  of  Xg  (elevator 
deflection)  and  chose  the  rudder  deflection  x^  =  0.  The  paths  Xg  >  Xg  =  0 
with  •>*  -0.0061771  have  two  limit  points,  for  >  Xg  >  Xg  =  0,  with 

aj2  -0.012498  a  third  limit  point  appears,  and  for  0)2  >  Xg,  Xg  =  0  only  one 
limit  point  remains.  A  similar  picture  arises  for  negative  roll  rates. 

In  all  cases  the  programs  easily  detected  and  computed  the  various  limit 

points  (see  Table  5).  But  the  example  also  shows  that  even  with  a  large  number  of  search 

paths  it  is  difficult  to  provide  a  full  picture  of  the  location  of  the  critical 

boundary,  that  is,  of  the  curves  of  limit  points  with  respect  to  x^  for  Xg  =  0 
and  varying  Xg.x^.  In  Figure  3  the  corresponding  branches  of  limit  point 
curves  are  shown  as  dotted  lines.  They  were  obtained  with  a  code  for  the 
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earlier  mentioned  task  (1.3)  (iv)  (see  [28]). 


"l 

^2 

^2 

"4 

"5 

"6 

^7 

1 

2.9649 

0.82557 

0.073661 

0.041309 

0.26735 

-0.05 

0.50481 

2 

1 

2.8174 

-0.17629 

0.089926 

0.026429 

-0.071476 

-0.008 

-0.20497 

3 

3.7579 

-0.65542 

0.38658 

0.092521 

-0.19867 

-0.008 

0.006208 

4 

4.1638 

0.089131 

0.094806 

0.022889 

0.016232 

-0.008 

-0.37766 

5 

2.5873 

[ 

-0.22355 

0.054682 

0.013676 

-0.091687 

-0.18691 

6  1 

3.9005 

-1.1482 

0.58156 

0.13352 

-0.32859 

0.51016 

1 

7 

2,2992 

-1.4102 

-0.061849 

-0.079009 

-0.58630 

0.1 

-0.68972 

8 

4.4565 

-4.4909 

1.6164 

0.33091 

-1 .0857 

0.1 

10.0212 

Table  5 

Example  3.  As  an  example  for  the  numerical  investigation  of  the  equilibrium 

surface  of  a  mechanical  structure,  we  consider  a  clamped,  thin,  shallow, 

circular  arch  which  has  been  used  as  a  test  case  by  various  authors  (see  eg. 

[12],  [16],  [36]),  Cet  U  and  W  be  the  radial  and  axial  displacements,  R 

the  arch  radius,  A  the  cross-sectional  area,  H  the  thickness,  and  E 

Young's  modulus.  With  the  dimensional  displacements  u  =  U/H,  w  »  W/H,  the 

2 

total  potential  energy  —  non-dimensional i zed  by  dividing  by  EAR(H/R)  — 
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is  given  by 


(6.2) 


{[(^  .  u)  +  1  -  + 


(4)'- 

de'^ 


0(2  p  u}de. 


Here  p  =  p(9)  is  the  dimensionless  radial  load,  and  a  i,  a2  are  dimensionless 
constants.  Each  end  is  assumed  to  be  pinned,  that  is,  we  have  the  boundary 
condi ti ons 


(6.3) 


u(+  Sg)  =  0, 


w(+  0^)  =  0, 


(±  V 


0. 


The  finite  element  approximation  introduced  in  [361  was  applied.  More 
specifically,  we  used  a  uniform  mesh  with  eight  elements,  0^  =  15  and  the 
constants  *  3.8716  x  10“^,  oig  *  1.65504  x  lo"^  corresponding  to  the 

data  in  [16].  Moreover,  the  following  load  function  p  =  p(y,v)  was  chosen 

fu(l  +  7v),  for  element  4 


(6.4) 


P(u.v)  =  < 


u(l  -  v) ,  otherwise 


corresponding  to  a  base  load  8  *  u(l-v)  and  an  excess  load  8uv  in  element 
4  such  that  the  average  load  is  always  u. 

Several  curves  on  the  equilibrium  surface  corresponding  to  constant 
values  of  u  or  V  were  computed.  Figure  4  shows  the  projection  of  these 
curves  into  the  (6,5)-plane  where  6  represents  the  radial  displacement  of  the 
center  point.  For  uniform  loads,  that  is,  v  =  0,  we  encounter  two  bifurcation 
points  on  the  primary  curve  which  are  connected  by  two  "buckling"  curves  that 
have  the  same  projection  in  the  (6,6)-plane. 
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THE  PITCON  CODE 


SUBROUTINE  F-TTCnN(NUBR»l  I^^.ITfXITfKSTEPf  IPC»H»  tRET»  TMOn.  TPVT.  PTCNOOOl 

1  HMAX»HHIN»HFACT.A8SERR..  *•  :RR.R«0RK»ISIZE1  PTrN0002 

C  PTCNOOO' 

C  PTCNOrtOfv 

C  PITCON. FOR  PTCNOOOa 

C  PTCNOOO? 

C  29  OCT  1981  VERSION  OF  PITCON»  PTCNOOOS 

C  THE  UNIVERSITY  OF  PITTSBURGH  CONTINUATION  PACKAGE.  PTCN0009 

C  THIS  VERSION  USES  SINGLE  PRECISION  AND  FUl.l.  NATRIX  STORAGE.  PTTNOOtO 

C  PTCNOOll 

C  PTC.N0012 

C  THIS  PACKAGE  NAS  PREPARED  WITH  THE  PARTIAL  SUPPORT  OF  PTCN0013 

C  THE  NATIONAL  SCIENCE  FOUNDATION*  UNDER  GRANT  HCS-78-0S299.  PTCN0014 

C  PTCN0015 

C  BY  MERNER  C.  RHEINBOLDT  AND  .JOHN  V.  BURKARDT  PTCN0016 

C  INSTITUTE  FOR  COHPUTATIONAL  HATHEHATICS  AND  APPLICATIONS  PTCNOOl? 

C  DEPARTMENT  OF  MATHEMATICS  AND  STATISTICS*  PTCN0013 

C  UNIVERSITY  OF  PITTSBURGH.  PITTSBURGH.  PA  152A1.  PTCN0019 

C  PTCN0020 

C  PTCNOOSl 

C  AN  EARLIER  VERSION  OF  THE  PACKAGE  NAS  URITTEN  IN  COOPERATION  WITH  PTCN0022 

C  PTCN002.i 

C  GEORGE  D  BYRNE*  PTCN0024 

C  COMPUTING  TECHNOLOGY  AND  SERVICES  PTCN0025 

C  EXXON  RESEARCH  AND  ENGINEERING  COMPANY  PTrNO026 

C  LINDEN.  MEM  JERSEY.  0703A  PTCN0027 

C  PTCN0028 

C  PTCN0029 

C  THIS  PACKAGE  COMPUTES  POINTS  ALONG  A  SOl.UTION  CURVE  OF  AN  PTCN0030 

C  UNDERDETERMINED  SYSTEM  OF  NONLINEAR  EQUATIONS  OF  THE  FORM  FX»0.  PTCN0031 

C  THE  CURVE  IS  SPECIFIED  TO  PASS  THROUGH  A  GIVEN  STARTING  SOLUTION  PTCN0032 

C  X  OF  THE  SYSTEM.  HERE  X  DENOTES  A  REAL  VECTOR  OF  NVAR  PTCN0033 

C  COMPONENTS  AND  FX  A  REAl.  VECTOR  OF  NVAR-l  COMPONENTS.  PTCNOO^X 

C  NORMALLY  EACH  CALL  TO  PITCON  PRODUCES  A  NEU  POINT  FURTHER  ALONG  PTCN0035 

C  THE  SOLUTION  CURVE  IN  A  USER-SPECIFIED  DIRECTION.  PTCN003e 

C  AN  OPTION  ALLOHS  THE  SEARCH  FOR  AND  COMPUTATION  OF  TARGET  POINTS.  PTCN0037 

C  THAT  IS.  SOLUTION  POINTS  X  FOR  MHICH  X<IT)  =  XIT  FOR  SOME  USER  PTCN0038 

C  SPECIFIED  VALUES  OF  IT  AND  XIT.  PTCN0039 

C  A  FURTHER  OPTION  ALLOWS  THE  SEARCH  FOR  AND  COMPUTATION  OF  LIMIT  PTCN0040 

C  POINTS  FOR  SPECIFIED  COORDINATE  LIM.  THAT  IS.  SOLUTION  POINTS  FOR  PTCN0041 

C  WHICH  THE  l.IM-TH  COMPONENT  OF  THE  TANGENT  VECTOR  IS  ZERO.  PTCNOO-i: 

C  PTrN0043 

C  PTCNf)044 

C  EXPLANATIONS  OF  THE  ALGORITHMS  USED  IN  THIS  PACKAGE  MAY  PTCN004f. 

C  BE  FOUND  IN  THE  FOLl.OHING  REFERENCES?  PTCN0046 

C  PTCN0047 

C  PTCNO048 

C  WFRNER  RHEINBOi  DT.  PTrN0049 

C  SOLUTION  FIELD 'of  NONLINEAR  EQUATIONS  AND  CONTINUATION  METHODS  PTCN0050 

C  SIAM  JOURNAL  OF  NUMERICAL  ANALYSIS.  17.  1980.  PP  221-237  PTCN0051 

C  PTCN0052 

C  COR  DEM  HEIJFR  AND  WERNER  RHEINBOl.DT.  PTCN0053 

C  ON  STEPLEN6TH  ALGORITHMS  FOR  A  CLASS  OF  CONTINUATION  METHODS.  PTCN0054 

C  SIAM  .JOURNAL  OF  NUMERICAL  ANAl.YSIS  18.  1981.  PP  925-947  PTCN0055 

C  PTCN005A 

C  WERMER  RHEINBOLDT.  PTCN0057 

C  NUMERICAL  ANAl.YSIS  OF  CONTINUATION  METHODS  FOR  NONl.INEAR  PTCN0058 

C  STRUCTURAL  PROBLEMS.  PTCN005? 

C  COMPUTERS  AND  STRUCTURES.  13.  1981.  PP  103-114  PTCNOO.SO 

C  PTCN00<S1 

C  PTCNOOA2 

C  OUTLINE  OF  THE  MATHEMATICAL  PROCEDURE  PTCN00A3 

C  PTCN00A4 

C  PTCNOO^S 

r  THE  FUNCTION  F  DEFINING  THE  SYSTEM  OF  EQUATIONS  IS  SUPPOSED  TO  PTCNOOSc 
C  BE  CONTINUOUSLY  niFFERENTIABI.E.  THE  REGULARITY  SET  R(F)  OF  F  IS  OTCNCOA? 

C  THF  SET  OF  ALL  POINTS  '<  IN  NUAR-DIMENSIONAl.  SPACE  WHERE  THE  pTCNOO.SB 

C  derivative  DF^X)  has  maximal  rank,  the  STARTING  POINT  IS  ASSUMED  PTCNOO.s? 

C  TO  BE  IN  R<FJ.  FOR  ANY  INDEX  IP  WITH  1 .LE.TP.LF.NVAR.  LET  PTCNLOTO 


OUTLINE  OF  THE  MATHEMATICAL  PROCEDURE 


THE  FUNCTION  F  DEFINING  THE  SYSTEM  OF  EQUATIONS  IS  SUPPOSED  TO 
BE  CONTINUOUSLY  niFFERENTIABI.E.  THE  REGULARITY  SET  R(F)  OF  F  IS 
THE  SET  OF  ALL  POINTS  '<  IN  NUAR-DIMENSIONAl.  SPACE  WHERE  THE 

derivative  dfix)  has  maximal  rank,  the  starting  point  is  assumed 

TO  BE  IN  R<FJ.  FOR  ANY  INDEX  IP  WITH  1 .LE.TP.LF.NVAR.  LET 


oooo.Toooooooooriooooonoooooooooooooooooorjrjoorjnooooonno  ooooonos-iooooooorjo 
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FACX.tp)  BE  THE  FUNCTION  OBTAIMEU  BY  AllGHENTING  F  JITH 
AN  NWAR-TH  SCALAR  FUNCTION  X<IP)-B  FOR  SOME^NUMBER  B.  FOR  ANY 
X  IN  R<F)  THERE  IS  AT  LEAST  ONE  INDEX  IP  SUCH  THAT  THE  DERIVATIVE 
DFAIXflP)  OF  FA  IS  NOMSINOULAR.  WITH  SUCH  AN  IP.  AND  A  niRECTTON 
DIR  (=+l  OR  -I).  THE  TANGENT  OF  FA.  TAN(X).  IS  UNTOUELY  DEFINED  BY: 

tan:*  SOJUTION  of  <  DFA<X.IP)*TAN=E(NVAR)  > 
tan:*  TAN/<EIJCLIDEAN  NORN  OF  TAN) 

SN:*  DIRBSIGN  (determinant  (DFA(X.IP)  )  ) 
tan:*  snbtan 

here  E(I)  is  the  I-TH  BASIS  VECTOR  IN  NVAR-SPACE. 

THE  PROCESS  USES  A  LOCAL  PARAMETERIZATION  OF  THE  CURVE. 

NORMALLY  THE  CONTINUATION  PARAMETER  INDEX  IP  IS  CHOSEN  AS 
THAT  INDEX  FOR  WHICH  ABS(TAN<X) ( IP) )  IS  MAXIMAL.  BUT  IN  THE 
CASE  OF  CERTAIN  CURVATURE  CHANGES  WHERE  IT  APPEARS  THAT 
A  LIMIT  POINT  FOR  THIS  CHOICE  FOR  IP  IS  APPROACHING. 

OTHER  CHfilCES  FOR  IP  MAY  BE  USED. 

PREDICTION  TAKES  PLACE  ALONG  THE  EULER  LINE  X+H*TAN.  THE 
STEPLENGTH  ALGORITHM  TAKES  INTO  ACCOUNT  THE  QUALITY  OF  THE 
CORRECTOR  ITERATION  AT  THE  LAST  POINT  AND  A  PREDICTION  OF  THE 
CHANGE  IN  CURVATURE.  THE  TANGENTIAL  STEPSIZE  USED  IN  PREDICTION  IS 
CHOSEN  SO  AS  TO  ACHIEVE  APPROXIMATELY  THE  PREDICTED  SECANT  STEPSIZE 
AFTER  CORRECTION  IS  DONE. 

THE  CORRECTOR  ITERATION  STARTS  FROM  THE  PREDICTED  POINT  AND  SOLVES 
THE  AUGMENTED  SYSTEM  FA(X.IP)=0  WITH  THE  VALUE  OF  THE  SCAl.AR 
B  EQUAL.  TO  THE  IP-TH  COMPONENT  OF  THE  PREDICTED  POINT.  THE  USER 
CAN  SPECIFY  AS  CORRECTOR  ITERATION  EITHER  A  FULL  NEWTON  PROCESS. 

OR  A  MODIFIED  NEWTON  PROCESS  WITH  FIXED  .JACOBIAN  DFA  EVALUATED  AT 
THE  PREDICTED  POINT. 


OUTLINE  OF  THE  COHPUTATIONAl.  ALGORITHM 


DURING  THE  FOLLOWING  DESCRIPTION.  WE  WILL  AS-SUME  THAT  WE 
HAVE  ENTERED  THE  CONTINLJATION  LOOP  WITH  AN  OLD  POINT  XL. 

A  CURRENT  POINT  XC.  THE  TANGENT  Tl.  AT  XL.  AND  CERTAIN  SCAl.AR 
miANTTTIES  ASSOCIATED  WITH  THESE  VECTORS.  WE  WILL  CHECK 
FIRST  FOR  ANY  TARGET  OR  LIMIT  POINTS  BETWEEN  XL  AND  XC. 

THEN  PROCEED  TO  COMPUTE  A  NEW  CONTINUATION  POINT  XF. 

THESE  NAMES  ARE  NOT  IN  PRECISE  ACCORDANCE  WITH  THE  STORAGE 
ARRANGEMENTS  UNTIL  THE  END  OF  A  CONTINLJATION  STEP. 


STEP  i: 


FOR  KSTEP.GT.O.  THE  CODE  GOES  TO  STEP  2. 

ON  THE  FIRST  CALL  TO  PITCON  FOR  A  GIVEN  PROBLEM  (KSTFP=-\ 

OR  KSTEP*0)  PROBLEM-DEPENDENT  CONSTANTS  ARE  SET 
AND  USER  CONTROL  PARAMETERS  ARE  LOADED  OR  DEFAULTS  USED. 

IF  (KSTEP.EO.O).  THE  PROGRAM  PROCEEDS  TO  STEP  2. 

IF  (KSTEP.EO.-L).  THE' USER  REQUESTS  THAT  THE  INPUT  STARTING  PTCN012.1 
POINT  XR  BE  CHECKED  FOR  THE  CONDITION  PTCN0122 

ABS<F(XR)).l^.(ABSERR/2).  IF  THIS  IS  NOT  THE  CASE.  NEWTON 'SPTCN012S 
METHOD  IS  APPLIED  TO  THE  POINT  XR  UNTIL  THE  ERROR  CONDITION  PTrN0t24 


PTCNOOTJ 

PTfNOOZZ 

PTCNOO?' 

•^TCNOO^-J 

PTrNO'ITf, 

PTrNo:*c 

PTCNOO/Z 

PTfUCOTS 

PTCN0079 

PTCNOuAO 

PTCN008) 

PTCN0032 

PTCNOOB3 

PTCN00fl4 

PTCNOOGS 

PTCN0086 

PTCN0087 

PTCN00S8 

PTCN0089 

PTC,N0090 

PTCN0091 

PTCN0092 

PTrN()093 

pTp(vrjrtP4 

PTCN0695 

pti:noo96 

PTCN0O97 

PTCN009S 

PTCN0099 

PTCNOIVO 

PTCNOLOI 

F'TCUC.tC2 

PTCN0I03 

PTCN0.104 

PTCNOJOS 

PTCNOlOft 

PTCN0107 

PTCNOtOS 

PTCN0109 

PTCNOllO 

PTCNOIJI 

PTCNO!12 

PTCN0n3 

PTCN0.H4 

PTCN0115 

PTCN01.L4 

PTCN0)17 

PTCNOllS 

PTCNOLl? 

PTCN0120 


IS  SATISFIED.  OR  A  FAILURE  HAS  OCCURRED.  AN  UNIMPROVABLE 
POINT  RESUI.TS  IN  A  RETURN  OF  IRET*-.4. 

IF  THE  STARTING  POINT  XR  WAS  IHPROVED.  THE  PROGRAM  RETURNS 
WITH  IRET*0  AND  KSTEP-0. 

IF  (KSTEP.EO.O).  THE  CONTINUATION  LOOP  BEGINS  WITH  THE 
STARTING  POINT  XR  STORED  IN  XL  AND  XC.  THE  STEPSIZE  HTANCF 
SET  TO  THE  INPUT  VALUE  OF  H.  AND  THE 
CONTINUATION  PARAMETER  SET  TO  THE  INPUT  VALUE  OF  IPC. 

FOR  KSTEP.GT.O.  THESE  OUANTITIES  ARE  COMPUTED  AND 
UPDATED  BY  THE  PROGRAM  ITSELF. 


STEP 


PTCN012S 
PTr.N0l26 
PTCNOtZ? 
PTCN0t2S 
PTCN0129 
PTCN0130 
PTCNOni 
PTr.N0l32 
PTCNOl  .-3 
PTCNOt34 
PTCN0(3S 
PTCNOl. 36 
PTCN0137 


TARGET  POINT  CHECK.  IF  (IT.NE.O).  A  TARGFT  POINT  IS 
DESIRED.  THE  VALUES  OF  XL (IT)  AND  XC(IT)  ARE  COMPARED  TO 
XIT.  IF  THE  TARGET  POINT  IS  BETWEEN  XI.  AND  XC.  THE  PROGRAM  PTCNOl 38 
COMPUTES  THE  TARGET  POINT.  SETS  IRET*1.  AND  PTCN013'? 

RETURNS.  TEMPORARILY  INTERRUPTING  NORMAL  CONTtNllATTON.  PTC.-JOJjO 


STEP  V, 


STEP  4: 


STEP  5J 


STEP  a: 


F'TrN0141 

tangent  and  l.OCAL  CONTINUATION  PARAHETER  CALCULATION.  IF  THEFTCNOMC 
LOOP  UAS  SUSPENDED  AT  THE  LAST  CALI.  TO  PITCON  TO  ALLOW  THE  PTCNOt43 
RETURN  OF  A  LIMIT  POINT>  THEN  THE  TANGENT  HAS  ALREADY  SEEN  PTCNOIAA 
CALCULATED  AND  A  LIMIT  POINT  CHECK  IS  SUPERFLUOUS^  SO  PTCNOIAE 
THE  PROGRAM  SKIPS  TO  STEP  5.  i 
OTHERUISEi  A  VECTOR  IN  THE  TANGENT  PLANE  AT  XC  IS  COMPUTED.  PTCNOJ-)? 
SUPPOSE  THAT  THE  PREVIOUS  CONTINUATION  PARAMETER  INDEX  PTCNOt43 
UAS  IPLi  WHERE  ON  THE  FIRST  STEP  IPl.  IS  USER-SUPPI. lED.  PTCNOIA? 
THE  NEW  TANGENT  IS  NORMALIZED^  AND  THE  IPL-TH  COMPONENT  PTCN0150 
IS  FORCED  TO  HAVE  THE  SAME  SIGN  AS  THE  IPL-TH  COMPONENT  PTCN0151 
OF  THE  PREVIOUS  TANGENT  (OR  ON  FIRST  STEP.  THE  SAME  ?TrN0152 
SIGN  AS  THE  USER  INPUT  DIRECTION  DIR.)  THEN  THE  LOCAL  PTCNOlfS 
CONTINUATION  PARAMETER  IPC  IS  DETERMINED.  IPC  IS  TO  rHEPTCN0134 
LOCATION  OF  THE  LARGEST  COMPONENT  OF  THE  TANGENT  VECTOR  PTCNOtSS 
UNLESS  A  LIMIT  POINT  FOR  THIS  CHOICE  APPEARS  TO  3E  PTCNOif!* 
APPROACHlNGf  IN  WHICH  CASE  THE  LOCATION  OF  THE  SECOND  PTCN0JS7 
LARGEST  COMPONENT  MAY  BE  TRIED.  PTCNOISS 
ONCE  IPC  IS  SET.  THE  RUANTITIES  TCIPC.  TCLIM.  HSFCLC.  ALPHLCPTCN0159 
AND  DIR  ARE  COMPUTED.  UHO.SE  MEANINGS  ARE  EXPLAINED  BELOW.  PTCN0160 

PTCNOlol 

LIMIT  POINT  CHECK.  IF  LIM.NE.O.  THE  LIM-TH  COMPONENTS  PTCNOtft2 
OF  THE  OLD  AND  NEW  TANGENTS  ARE  COMPARED.  IF  THESE  DIFFER  PTCNOlcS.T 
IN  SIGN.  A  LIMIT  POINT  LIES  BETWEEN  XL  AND  XC.  THE  PROGRAM  prCMl64 
ATTEMPTS  TO  FIND  THIS  LIMIT  POINT.  IF  FOUND.  IT  STORES  PTCNOIaS 
THE  LIMIT  POINT  IN  XR.  THE  TANGENT  AT  XR  IN  TL.  SETS  IRET=2.PTCN016o 


AND  RETURNS.  TEMPORARILY  INTERRUPTING  THE  NORMAL  LOOP. 


THE 


STEP  7: 


STEP  LENGTH  COMPUTATION.  THE  PROGRAM  COMPUTES  HTANCF 
STEPSIZE  TO  BE  USED  ALONG  THE  TANGENT  TO  OBTAIN  THE 
PREDICTED  POINT  XPREDJsXCFHTANCF»TC..  THE  STARTING  POINT 
FOR  THE  CORRECTOR  PROCESS.  IN  COMPUTING  HTANCF.  CERTAIN 
CURVATURE  AND  STEPSIZE  DATA  ARE  UPDATED. 

PPEmCTION  AND  CORRECTION  STEP.  WITH  THE  PREDICTED  POINT 
XPRED*XCFHfANCFtTC  AS  A  STARTING  POINT.  THE  CORRECTOR 
PROCESS  IS  APPLIED  TO  CORRECT  THE  POINT  UNTIL 
ABS<F(XCOR)).l.E.ABSERR  AMD  XSTEP.  THE  LAST  CORRECTOR  STEP. 
SATISFIES  XSTEP. LE.ABSERR+REl.ERR*ABS(XCOR) . 

IF  THE  SIZE  OF  A  CORRECTOR  STEP  IS  TOO  LARGE. 

OR  IF  A  CORRECTION  STEP  INCREASES  THE  FUNCTION  VALUE.  OR 
THE  MAXIMUM  NUMBER  OF  STEPS  ARE  TAKEN  WITHOUT  CONVERGENCE. 
THE  STEPSIZE  HTANCF  IS  REDUCED  AND  THE  CORRECTOR  STEP  IS 
ATTEMPTED  AGAIN.  IF  THE  STEPSIZE  SHRINKS  BEl.OW  HMIN.  THE 
PROGRAM  SETS  AN  ERROR  FLAG  AND  RETURNS. 

STORING  INFORMATION  BEFORE  RETURN.  AFTER  A  SUCCESSFUL 
CONTINUATION  STEP.  THE  PROGRAM  REARRANGES  ITS  STORAGE  SO 


PTCN01<f.7 

PTCNOIXS 

PTCN0169 

PTCN0170 

PTCN0171 

PTrN0172 

PTCN0173 

PTCN0(74 

PTCN0I75 

PTCN0t76 

PTrN0177 

PTCN0179 

PTCN0t79 

PTCNOIRO 

PTCN0181 

PTCN0182 

PTCN0IB3 

PTCN0134 

PTCNOISS 

PTCNOIS* 

PTCN01S7 

PTr.N01.33 


THAT  THE  ENTRIES  CORRESPONDING  TO  XC  AND  XF  HOLD  THE  PROPER  PTCNOtS? 


DATA.  COMPUTES  CORDIS.  THE  SIZE  OF  THE  CORRECTION  TO  THE 
PREDICTED  POINT.  AND  MODIFIES  CORDIS  TO  A  VALUE  THAT  WOULD 
CORRESPOND  TO  AN  OPTIMAL  NUMBER  OF  CORRECTOR  STEPS. 


ON  NORMAL  RETURN,  THE  VECTOR  XR  (THE  FIRST  NVAR  ENTRIES  OF  RUORK>. 
CONTAINS  A  SOLUTION  POINT  ON  THE  CURVE  F(XR)=0,  AND  IS  EITHER 
A  CONTINUATION  POINT.  A  TARGET  POINT.  OR  A  LIMIT  POINT.  WHICH 
IS  INDICATED  BY  THE  VAL.UE  OF  IRET. 

IF  IRET  IS  NEGATIVE,  AN  ERROR  HAS  OCCURRED.  IF  A  LIMIT  POINT  IS 
RETURNED,  THE  TANGENT  VECTOR  AT  THE  LIMIT  POINT  IS  CONTAINED  IN  THE 
LOCATION  TL  IN  RWORK.  ON  FIRST  CALL.  THE  USER  MUST  SET  SOME  OF  THE 
SCAI.AR  PARAMETERS,  AND  THE  STARTING  POINT  XR.  THEREAFTER,  ONLY  IT 
AND  XIT  SHOULD  BE  CHANGED  BY  THE  USER  DURING  A  PROBLEM  RUN. 

IF  A  NEW  PROBLEM  IS  TO  BE  RUN  (WHETHER  A  DIFFERENT  FUNCTION.  OR 
SAME  FUNCTION  WITH  DIFFERENT  STARTING  POINT  OR  ERROR  CONTROLS) 

THE  PROGRAM  MAY  BE  RESET  BY  USING  KSTEP=-I  OR  0,  AT  WHICH  TlnE  THE 
SCALARS  AND  THE  POINT  XR  MUST  BE  SET  AGAIN.  NOTE  THAT  IN  THIS  CASE 
THE  STATISTICAL  DATA  IN  THE  COMMON  BLOCKS  /COUNTl/  AND  /COUNT:. 

WILL  BE  RESET  TO  0  AS  HELL. 


PTCN0190 
PTCN0191 
PTCN0192 
PTCN0193 
PrCN0t?4 
PTCN019S 
PTCNOl?* 
PTCN0t97 
PTCNOl'78 
PTCN0199 
PTCN0200 
PTCN020I 
PTCN02D2 
PTCN0203 
THEPTCNOTOJ 
FTCN0205 
PTCNOT''* 
P’’CN0207 
PTC  NO'  -.3 
PTC NOTmc 
^Tj-NOri', 
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nFFINITIONS  AND  nEFAIU.TS  OF  PTTCOlV  PARArtETFRS 


MVAR  =  THE  NUMBER  OF  VARIABLES  IN  THE  NONLINEAR  SYSTEM.  NVAR  IS 

THE  DIMENSION  OF  THE  PIVOT  VECTOR  IPVTf  AND  THE  SI7E  OF  THE 
VECTORS  XRi  XC.  XF*  Tl.  AND  TC  MHICH  ARE  CONTAINED  IN  RUORK . 

RUORK  ALSO  CONTAINS  STORAGE  FOR  THE  MATRIX  FPRYN  WHICH 
IS  OF  SIZE  NVAR  X  NVAR. 

NVAR  MUST  BE  GREATER  THAN  1.  AND  MUST  WOT  BE  CHANGED  DURING 
THE  COURSE  OF  A  PROBLEM  RUN.  NVAR  HAS  NO  DF.FAUl  T  VALUE. 

LIM  =  LIMIT  POINT  FLAG  AND  INDEX.  IF  (LIM.EO.O).  LIMIT  POINTS 
ARE  NOT  TO  BE  LOOKED  FOR.  OTHERWISE f  THE  USER  SHOULD  SET 
l.IM  TO  THE  INDEX  OF  THE  VARIABLE  FOR  WHICH  L  IMIT 
POINTS  ARE  TO  BE  SOUGHT.  LIM  DEFAULTS  TO  ZERO. 

LIM  MLIST  SATISFY  O.l.E.LIM.LE.NVAR. 

IT  =  TARGET  POINT  FLAG  AND  INDEX.  IF  <IT.EQ.0)f  TARGET  POINTS 
ARE  NOT  TO  BE  L.flOKED  FOR.  OTHERWISE,  THE  USER  SHOULD  SET 
IT  TO  THE  INDEX  OF  THE  VARIABLE  FOR  WHICH  TARGET 
VAI.UES  XIT  ARE  DESIRED.  IT  HAS  THE  DEFALil.T  VALUE  ZERO. 

IT  HAY  DE  RESET  BY  THE  USER  AT  ANY  TINE  DURING  A  RUN. 

IT  MUST  SATISFY  O.LE.IT.LE.NVAR. 

XIT  »  THE  VAl.UE  OF  THE  TARGET  VECTOR  COMPONENT  SOUGHT,  IF  IT.NE.O.  PTCN02S3 
TARGET  POINTS  XR  SATISFY  XR(IT)*XIT.  XIT  HAS  NO  DEFAULT.  PTCNOZZA 

KSTEP  »  THE  NUMBER  OF  CONTINLIATION  STEPS  TAKEN.  THIS  DOES  NOT  FTCN0Z3S 

INCl.UDE  FAILURES,  TARGET  POINTS  OR  LIMIT  POINTS.  THE  PR0GRAMPTCN0236 
INCRENENTS  KSTEP  EACH  TIME  A  NEW  POINT  XE  IS  COMPUTED.  PTCN0237 

ON  THE  FIRST  CALL  TO  PITCON  FOR  A  PARTICULAR  PROBl.EM,  THE  PTCN0238 

USER  SHOULD  SET  KSTEP  TO  0  OR  -I.  IE  KSTEP=-1,  THE  PROGRAM  PTCN0239 
HILL  CHECK  THE  ACCURACY  OF  THE  STARTING  POINT  XR,  AND  IF  PTCN0240 
NECESSARY,  ATTEMPT  TO  CORRECT  IT  USING  NEWTON'S  METHOD.  FTCN0241 

IF  KSTEP=0,  THE  PROGRAM  PERFORMS  NO  CHECK  ON  THE  START 'NG  PTCN0242 
POINT,  AND  PROCEEDS  TO  THE  CONTINUATION  LOOP.  IF  THE  USER  PTCN0243 
WISHES  TO  RUN  A  DIFFERENT  PROBLEM  THEN  A  CALL.  TO  PITCON  PTCN0244 

WITH  KSTEP»-1  OR  0  WILL  RESET  THE  CODE,  DESTROYING  THE  PTCN0245 

INFORMATION  FROM  THE  PREVIOLIS  RUN.  KSTEP  DEFAULTS  TO  -t.  PTCN0246 

IPC  =  THE  COMPONENT  OF  XC  TO  BE  USED  AS  CONTINUATION  PARAMETER.  PTCN0247 

ON  THE  FIRST  CALL  ONLY,  THE  USER  OUGHT  TO  SET  IPC  OR  ALLOW  PTCN0748 

THE  DEFAULT  VALUE  IPC=NVAR.  AFTER  THE  FIRST  CALL,  THE  PTCN024'’ 

DETERMINATION  OF  IPC  IS  DONE  BY  THE  PROGRAM  USING  INFORMATIONFTCNO^SO 


r  rr;jO;i  : 

PTCNOZIZ 
FTrN0?.l4 
PTr,N0215 
PTCNOZli 
PTCNOn  7 
PTCNOZIS 
PTCN073'? 
PTrNr)220 
PTCNOZn 
pTCNo::: 
PTCN0223 
PTCNOrZ^ 
FTCN0725 
PTCNOZZc 
PTr.N0r7.7 
PrCN0223 
PTCNOJ'? 
PTCNGZZO 
PTCN023] 
FTCNO'32 


ABOUT  THE  TANGENT  VECTOR  AT  XC. 

H  =  SUOGESTED  STARTING  STEP  SIZE  Al.ONG  THE  TANGENT  TO  THE  CURVE. 
IF  H=0.0  ON  THE  INITIAL  CALL,  H  DEFAULTS  TO  (HMAX+HMIN) /2 
IF  H  IS  NEGATIVE  ON  THE  FIRST  CAI.L,  THE  MINUS  SIGN 
IS  ABSORBED  BY  DIR  AND  INDICATES  THAT  THE  DIRECTION  OF  THE 
FIRST  STEP  SHOULD  BE  IN  THE  NEGATIVE  IPC  DIRECTION. 

AFTER  THE  FIRST  STEP,  STEPSIZE  IS  CONTROLLED  BY  THE  PROGRAM. 
UPON  RETURN  WITH  A  CONTINUATION  POINT  (IRET=0),  H  IS 
OVERWRITTEN  BY  HTANCF,  THE  STEPSIZE  USED  IN  REACHING  THE 
NEW  POINT, 

IRET  »  A  RETURN  FLAG  TO  INDICATE  ERRORS  OR  THE  TYPE  OF  POINT 
RETURNED  IN  XR.  NONNEGATIVE  VAI.UES  OF  IRET  REPRESENT 
NORMAL  RETURNS.  HERBTIVJi  VAJ-UES  OF  IRET  INDICATE  THAT 
SOME  ERROR  OR  DIFFICUL.TY  HAS  BEEN  ENCOLJNTERED .  VAI.UES 
OF  IRET  BETWEEN  -1  AND  -4  ARE  SIMPLY  REPORTS  THAT  AN 
ATTEMPT  TO  COMPUTE  A  LIMIT  OR  TARGET  POINT  FAILED.  THERE 
DO  NOT  AFFECT  FURTHER  CONTINUATION  STEPS,  AND  THE  USER  NEED 
NOT  MODIFY  ANY  VARIABLES  BEFORE  PROCEEDING. 

VALUES  OF  IRET  OF  -5  AND  -A  REFER  TO  DANGEROUS  SITUATIONS 
THAT  MAY  BE  CORRECTABLE. 

VALUES  OF  IRET  FROM  -7  TO  -10  ARE  SERIOUS,  FATAL  ERRORS. 

THE  USER  SHOLJLD  HALT  THE  PROGRAM  AND  EXAMINE  HIS  INPUT 
AND  THE  INTERIM  RESULTS, 

IRET  SHOULD  BE  ZERO  ON  THE  FIRST  CALL.  FOR  A  PROBl.EM. 

THE  SPECIFIC  VALUES  OF  IRET  AND  THEIR  MEANINGS  ARE: 


IRETx2; 


NORMAL  RETURN  WITH  LIMIT  POINT  IN  XR  AND  TANGENT 
AT  XR  CONTAINED  IN  TL. 


IRFTat:  NORMAL  RETURN  WITH  TARGET  POINT  IN  XR. 


PTCWOZSJ 

PTCN02^: 

PTCN02S3 

PTCN02S4 

PTCN02S5 

PTCNOZSf 

PTrN02S7 

PTCW0258 

PTCN02S9 

PTCWOIoO 

PTCN02A1 

FTCN0762 

PTCN02A3 

PTr.N0264 

PTr,N02A5 

PTCN0266 

PTCNOZA? 

P*rN0268 

PTCN02A9 

PTCWOZZO 

PTCN027I 

PTrN0772 

PTCN0273 

PTrNDr74 

PTCN027S 

?TrN027a 

PTCN027'’ 

PTCN0278 

PTCW027« 

:rrwo:gc 
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IRFTsO:  NORMAL  RETURN  WITH  NEW  CONTINUATION  POINT  IN  XR. 


IRET=-.t:  AN  ERROR  OCCURRED  DURING  COMPUTATION  OF 


PTCNOrSl 
PTCN0?S2 
PTCNOrR? 
LIMIT  pnTNTPTrN0284 
PTCNO.tPS 


IRET=-.?: 


iret=-t; 


IRET=-4: 


IRET=-S: 


CORECT  CALLED  FOR  TARGET  POINT  CALCULATION  FAILED 
AFTER  KNMAX  ITERATIONS. 


SOLVE  WAS  CALLED  DY  CORECT  FOR  TARGET  POINT 
CALCULATION^  AND  FAILED.  (MATRIX  ELIMINATION 
ZERO  PIVOT). 


UNACCEPTABLE  CORRECTOR 
CALCULATION. 


STEP  IN  TARGET  POINT 


PREDICTION  STEP  HTANCF  IS  t.ESS 
BECAUSE  OF  REPEATED  FAILURE  OF 
CONSEQUENT  STEPSI7E  REDUCTION. 

HMIN.  OR  SWITCH  FROM  IM0D=1  TO 
ABSERR  AND  REl.ERR.  BUT  »E  AWARE  THAT  REPEATED 
STEPSI2E  REDUCTIONS  MAY  INDICATE  AN  INTRACTAPI.E 
FUNCTION. 


THAN  HMIN,  PERHAPS 
CORECT,  AND 
USER  MIGHT  REDUCE 
IM0D=O,  OR  INCREASE 


iret=-a: 


FUNCTION  VALUE  FNRMXF  OF  INPUT  XR  IS 
AND  COUl.D  MOT  BE  IMPROVED  BY  CORECT. 


TOO  LARGE 
USER 


PTCNnrS6 
PTCNOrP? 
PTC.N0:’S8 
PTCN02P9 
F0UNDPTCN0290 
PTCN029I 
PTCN0292 
PTCN029S 
PTCN0294 
PTCN0295 
PTCN029& 
PTCN0297 
PTCN0298 
PTCN0299 
PTCNOSOC 
PTCNOSOt 
PTC.N0302 
PTCNOSOS 
PTCN0.TA4 
PTCN0S05 


MIGHT  RECOVER  BY  RELAXING  ERROR  CONTROLS,  IMPR0VINGPTCN0S06 


STARTING  POINT  XR,  OR  CHANGING  VALUE  OF  IPC. 


IRET=-7:  SOLVE  FAILED  IN  A  CALL  FROM  TANGNT. 

IRET=-8:  SOLVE  FAILED  IN  A  CALL  FROM  CORECT. 

IRET=-9{  THE  TANGENT  VECTOR  TC  AT  XC  IS  ZERO. 

IRET=-10J  IMPROPER  INPUT,  NVAR.lE.t,  OR 
ISlZE.LT.<NVAR)«(NVAR+5),  OR 
PROGRAM  HAS  BEEN  CALLED  AGAIN  AFTER  FATAL  ERROR. 

IMOD  =  METHOD  FLAG  FOR  CORRECTOR  STEP,  SPECIFYING  TYPE  OF 
NEWTON  METHOD  TO  BE  USED. 

IM0D=0:  UPDATE  JACOBIAN  FOR  TANGENT  CALCULATION, 

UPDATE  JACOBIAN  FOR  EACH  CORRECTOR  STEP, 

IMOD*i:  UPDATE  JACOBIAN  FOR  TANGENT  CALCULATION, 

EVALUATE  JACOBIAN  AT  FIRST  CORRECTOR  STEP  ONLY. 


PTCN0307 
P  rCNOTOS 
PTCN0309 
PTCNOTIC 
PTCN0311 
PTCN0312 
PTCN0313 
PTCN0314 
PTCN03.t5 
PTCN03J6 
PTCN0317 
PTCNOSIS 
PTCN03.t9 
PTCN0320 
PTCNO.'ZI 
PTCN0322 
PTCN0323 
PTCN0324 
PTCN0325 
PTCN0326 
PTCN0327 
PTrN03?8 
PTCNO'29 
PTCN0330 
PTCN0331 
PTrN0332 
PTCN0333 
PTCN0334 
PTrN033S 
PTCN0336 
PTCN0337 


IPVT  =  INTEGER  VECTOR  USER  DECLARED  TO  BE  OF  SIZE  NUAR. 

USED  DURING  THE  MATRIX  FACTORIZATION  TO  STORE  PIVOT 
INFORMATION. 

HMAX  s  THE  MAXIMUM  STEP  SIZE.  IF  HMAX.IE.O.O  ON  INITIAL  CALL, 

HMAX  DEFAUl.TS  TO  SQRT(NVAR). 

HMIN  =  THE  MINIMUM  STEP  SIZE.  IF  HMIM.LE .SQRKEPMACH)  ON  INITIAL 
CALL,  HNIN  DEFAULTS  TO  SQRKEPMACH),  WHERE  EPMACH  IS  THE 
MACHINE  PRECISION  CONSTANT. 

HFACT  =  LIMIT  ON  STEPSIZE  CHANGE.  HSECLC  IS  THE  SECANT  STEPSIZE 
OF  THE  LAST  STEP,  AND  HTANCF  IS  THE  STEPSIZE  TO  BE  USED 
IN  OBTAINING  THE  PREDICTED  POINT,  THE  FOLl.OWING  RELATI0NSHIPPTCN0338 
MUST  BE  satisfied;  <HSFCLC/HFACT).LE, HTANCF. LE.<HSECLC«HEACT)PTCN0339 
IF  THE  CORRECTOR  STEP  FAILS,  THEN  HFACT  IS  ALSO  USED  TO  FTCN0340 
REDUCE  THE  PREDICTOR  STEP  HTANCF  TO  (HTANCF/HFACT)  PTCM0341 

IF  HFACT. LE, 1.0  ON  INITIAL  CALL,  HFACT  DEFAULTS  TO  3.0.  P’'CN0342 

ABSERR*  ABSOLUTE  ERROR  CONTROL.  IF  ABSERR=0.0  ON  INITIAL  CALL  PTCN0343 

ABSERR  DEFAULTS  TO  SORT(EPMACH)  FTrN0344 

RELERR*  RELATIVE  ERROR  CONTROL.  IF  RELERR=0.0  ON  INITIAL  CALL  PTCN0345 

RELERR  DEFAULTS  TO  SQRT(EPMACH)  PTrN034o 

RWORK  =  USER  DECLARED  VECTOR  OF  SIZE  ISIZE=NVAR*(NVAR+S) .  PTCN034? 

RWORK  STORES  FIVE  VECTORS  AND  AN  ARRAY  IN  THE  ORDER  PTCN03J8 

(XR,XC,XF,TI.,TC,FPRYM).  THEIR  BEGINNING  LOCATIONS  PTCNOSA? 

ARE  IXR*t,  IXCsNVAR+l,  IXFs.ZtNVARLl ,  ITL=3*NVAR+J .  PTCNOSSO 

ITC=4)(iNVAR+I,  IFP=S*NVAR+1  .  FPRYM  IS  AN  ARRAY  OF  ptCNOSSI 
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SIZE  MVAR  .X  MVAR.  THE  HEANINGS  OF  THESE  COMPONENTS  OE  PTCNOZ-'iC 

RUORK  ARE  DESCRIBED  BEl.OW.  THE  USER  SHOUI.D  SET  A  VAIHE  TO  XRPTCN0S53 
ON  FIRST  CALl.f  BUT  NO  OTHER  PORTIONS  OF  RUORK  SHOIJI  D  BE  PT:N0~5-» 

SET.  AFTER  THE  FIRST  CALL  FOR  A  PROBl.EN»  NO  ENTRIES  OF  RWORLFTCNOSSS 
SHOULD  BE  ALTERED.  PTCNOSSo 

(XR)  *  ON  FIRST  CALLf  A  USER  SUPPLIED  STARTING  POINT.  WHICH  NAY  BE  PTCNOSS? 

IMPROVED  BY  THE  PROGRAM  IF  KSTEP=-1,  ON  NORMAL  RETURN  PROn  PTCNOSSS 
PITCON.  XR  MILL  HOLD  THE  MOST  RECENTL  Y  FOUND  POINT.  WHETHER  PTCNO;??:? 
A  CONTINUATION  POINT.  TARGET  POINT.  OR  LIMIT  POINT.  PTCNPSoO 

<XC)  =  THE  PREVIOUS  CONTINUATION  POINT.  PTCNOSAt 

(XF)  =  THE  CURRENT  CONTINUATION  POINT.  P*CN07a2 

(TL)  =  PREVIOUS  VALUE  OF  TANGENT  VECTOR.  NOTE  THAT  THIS  CORRESPONDSPTCNO.lo.S 

TO  A  POINT  XL  WHICH  HAS  BEEN  DISCARDED.  ON  A  LIMIT  POINT  PTCN0364 
RETURN.  TL  UILL  CONTAIN  INSTEAD  THE  TANGENT  AT  THE  PTCN03o5 

l.IMIT  POINT.  ON  A  TARGET  POINT  RETURN.  Tl.  UILL  HAVE  PTCNO'-So 

BEEN  OVERWRITTEN  BY  THE  FUNCTION  VALUE  AT  THE  TARGET  POINT.  PTCN03<S7 
(TO  =  THE  TANr.ENT  VECTOR  AT  THE  PREVIOUS  CONTINUATION  POINT.  BTCN03-S8 

(FPRYN)=  MATRIX  STORAGE  AREA  FOR  SETTING  UP  AND  SOI  VI NR  THE  PTCN03«° 

LINEAR  SYSTEMS  INVOLVING  DFAIX.IP).  PTi:N0370 

ISIZE  »  USER  SET  DIMENSION  FOR  VECTOR  RWORK.  WHICH  MUST  BE  AT  PTCN037t 

LEAST  OF  SIZE  NVAR*<NVAR+5) .  PTrN0372 

PTCN0373 

PTCN0374 

NOMENCLATURE  FOR  STEP  DEPENDENT  VARIABLES  PTCN037S 


THE  PROGRAM  ACCUMULATES  INFORMATION  THAT  IS  ASSOCIATED  WITH 
SEVERAl.  PREVIOUS  CONTINUATION  POINTS  OR  THE  STEPS  MADE  BETWEEN 
THEM.  IN  INTERPRETING  THE  CODE  OR  ITS  OUTPUT.  IT  IS  IMPORTANT 
TO  KNOW  WHERE  SUCH  GUANTITIES  APPLY.  THE  FOLLOWING  DESCRIPTION 
OF  SOME  OF  THE  VARIABLES  IS  VALID  ONLY  LJPON  A  NORMAL  RETURN  WITH 
A  CONTINUATION  POINT. 

THE  POINTS  'XLL'  AND  'XL'  WILL.  HAVE  BEEN  DISCARDED  BY  THE  PROGRAM, 
BUT  SOME  OUANTITIES  ASSOCIATED  WITH  THEM  STILL  SURVIVE. 

OLIANTITIES  ASSOCIATED  WITH  STEP  FROM  'XLL'  TO  'XL': 

HSECLl.  »  SIZE  OF  SECANT  ERON  'XLL'  TO  'XL'.  EUCLIDEAN  NORM(XLL-Xl.) 

(JUANTITIES  ASSOCIATED  WITH  THE  POINT  'XL't 

IPL  «  THE  LOCATION  OF  THE  FIRST  OR  SECOND  LARGEST  COMPONENT 
OF  THE  TANGENT  VECTOR  AT  'XL'. 

TLLIM  *  VALUE  OF  LIM-TH  COMPONENT  OF  TANGENT  VECTOR  AT  'XL'. 

TL  =  TANGENT  VECTOR  AT  'XL'.  ALTHOUGH  LIMIT  OR  TARGET  POINT 
CALCUl.ATIONS  COULD  HAVE  OVERWRITTEN  THIS  VECTOR. 

QUANTITIES  ASSOCIATED  WITH  INTERVAL  FROM  'XL'  TO  XC; 

ALPHLC  s  THE  ANGLE  BETWEEN  THE  TANGENTS  AT  'XL '  AND  XC. 

CURVl.C  =  ESTIMATED  CJJRVATURE  BETWEEN  'XL'  AND  XC. 

H8ECI.C  »  SIZE  OF  SECANT  BETWEEN  'XL'  AND  XC.  EUCLIDEAN  NORM(Xl.-XC) 

QUANTITIES  ASSOCIATED  WITH  THE  POINT  XC: 

DETA  =  BINARY  MANTISSA  OF  DETERMINANT  OF  DFA(XC.IPL),  DIVIDED 
BY  IPL-TH  COMPONENT  OF  TANGENT  AT  XC. 

DIR  a  SIGN  OF  DETA,  DF-TERMINES  SENSE  OF  CONTINUATION. 
lEXP  a  BINARY  EXPONENT  OF  DETERMINANT  OF  DFA(XC, IPl.) . 

IPC  a  LOCATION  OF  FIRST  OR  SECOND  LARGFST  COMPONENT  OF  TANGENT 
VECTOR  AT  XC, 

TC  a  tangent  VECTOR  AT  XC. 

TCLIN  a  VAI.UE  OF  LIM-TH  COMPONENT  OF  TANGENT  AT  XC. 

TCIPC  a  value  OF  TC(IPC) 


PTCN0376 

PTCN0377 

PTCN037S 

PTCN0379 

PTCM0380 

PTCN03B1 

PTCN0382 

PTCN03B3 

PTCN03S4 

PTCN0385 

PTCN03R6 

PTCN0387 

PTCN0388 

PTCN03R9 

PTCN0390 

PTCN0391 

PTCN0392 

PTCN0393 

PTCN0394 

PTCN039S 

PTr.N0396 

PTCN0397 

PTrN039S 

FTr-N039» 

PTCN0400 

FTCN0401 

PrCN0402 

PTCN0403 

PTrN04O4 

PTCN0405 

PTCN0406 

PTCN0407 

PTCN0408 

PTr,N0409 

PTCN0410 

PTCN0411 

PTCN0412 

PTCN04t3 

PTCN041J 

PTCN04tS 


QUANTITIES  ASSOCIATED  WITH  THE  INTERVAL  FROM  XC  TO  XF*. 


PTr.N04U 

PTCN04I7 


PTCN0413 

CURVCF  a  estimated  CURVATURE  BETWEEN  XC  AND  XF.  PTCN04I9 

HTANCF  a  STEPSIZE  USED  ALONG  TANGENT  TO  GET  PREDICTED  POINT  PTCNO470 
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yHICH  UAS  i;URRFCTEn  TO  SOLUTION  POINT  XF. 

OUANTITIES  ASSOCIATED  WITH  THE  POINT  XFt 

CORDIS  *  SIZE  OF  THE  TOTAl.  CORRECTION  FROM  PREDICTED  POINT 
X=XCfHTANCF»TC  TO  CORRECTED  POINT  XF. 

NOTE  THAT  THIS  HAS  BEEN  MODIFIED  TO  AN  'OPTIMAL'  VALUE. 
CURVXF  =  A  PREDICTED  VALUE  OF  THE  CURVATURE  AT  XF . 

FNRMXF  =  MAXIMUM  NORM  OF  FUNCTION  VALUE  AT  XF . 

FPRYM  =  DFA<XFjIPC)  HAS  ACTUALLY  BEEN  LAST  EVALUATED  AT  THE 

PENULTIMATE  CORRECTOR  ITERATE  (IF  IMOD.NE.l).  IT  Will.  BE 
EVALUATED  AT  XF  AS  SOON  AS  THE  NEXT  LOOP  BEGINS  AND  THE 

XSTEP  =  SIZE^OF  THE  LASf  CORRECTOR  STEP  TAKEN  IN  CONVEROINO  TO  XF - 


SUBROUTINES  IN  THIS  PACKAGE 


PTCNOA?! 

PTCN0422 

PTCN042S 

PTCN0424 

PTCN042S 

PTCN042(i 

PTr.N0427 

PTrN0428 

PTrN0429 

PTCN0430 

PTCN04S1 

PTCN0432 

PTCN0437 

PTCN0434 

PTr.N0435 

FTCN043A 

PTCN0437 

PTCN0438 

PTCNOAS? 


PITCON<NVARjLIM»IT>XIT..KSTEP.IPC»H»IREY».IMOD»IPVT»  PTCN0440 

HMAX>HMIN>HFACTfABSERR>RELERRfRU0RKfISI7E)  PTCN0441 

PTCN0442 

DRIVING  ROUTINE  OF  CONTINUATION  CODE.  INITIALIZES  INFORMATION^  PTCN044.1 
DETERMINES  WHETHER  LIMIT»  TARGET  OR  CONTINUATION  POINT  HILL  PTCN0444 

BE  SOUGHT  THIS  STEP.  COMPUTES  STEPLENGTHS.  CONTROLS  CORRECTOR  PTCN0445 

PROCESS.  AND  HANDLES  ERROR  RETURNS.  PTCN0446 

PTCN0447 

CORECTXNVAR.X. IHOLD. WORK. lERR. IMOD> FPRYM. IPVT. ABSERR.RELERR.  PTCM0448 

XSTEP. NEQN.FNRM)  PTCN0449 

PTCN0450 

USFS  A  FORM  OF  NEWTON'S  METHOD  TO  SOLVE  THE  AUGMENTED  NONLINEAR  PTCM0451 
SYSTEM  FA(X)aO  WITH  AUGMENTING  EQUATION  X(IHOLD)=B.  THAT  IS.  X( IHOLD )PTCN0452 
IS  HELD  FIXED  DURING  THE  CORRECTION  PROCESS.  PTCN04S3 


PTCN04S4 

TANGNT  <  NVAR . XC . IPC . TC . I RET » I CAL  1. 1 FPR YM . IPVT . NEON . DET A . I EXP )  PTCN04S5 

PTCN0456 

applies  Al.GORITHM  DESCRIBED  ABOVE  TO  SOLVE  DFA(XC. IPL)*TC»E(NVAR)  PTCN0457 
AND  THEN  NORMALIZES  TANGENT  VECTOR.  CORRECTS  SIGN.  AND  SETS  PTCM0458 

IPC  AND  DIR.  PTCN04S9 

PTrN04A0 

ROOT  <  A . F A  f  B . FB . C . FC . KOUNT . IFL AG )  PTCN04 A I 

PTCN04o2 

ROOT  FINDER  USED  TO  LOCATE  LIMIT  POINT.  THIS  ROUTINE  IS  A  MODIFIED  PTCN04A3 
VERSION  OF  THE  FORTRAN  FUNCTION  ZERO  GIVEN  IN  THE  BOOK!  PTCN0464 

'ALGORITHMS  FOR  MINIMIZATION  WITHOUT  DERIVATIVES  PTCM04A5 

BY  RICHARD  P  BRENT.  PRENTICE  HAl.L.  1973.  PTCNOAAA 

PTCN04A7 

SOLVE  <  NVAR.X. Y . IP . DETA . lEXP  > lERR. ICALL . IMOD . FPRYM . IPVT )  PTCN04A8 

PTCN04A9 

SETS  UP  AND  SOLVES  THE  SYSTEM  DFA(X.IP)*Y(OUTPUT)=y(INPUT)  PTCN0470 

WHERE  DFA(X.IP)  IS  THE  JACOBIAN  OF  FA  AT  X.  PTCN0471 

AND  Y  IS  A  RIGHT  HAND  SIDE  SUPPLIED  BY  THE  CALLING  ROUTINE.  PTCN0472 


PTCM0473 

«*NOTE|[*  SUBROUTINE  SOl.VE  USES  FULL  MATRIX  STORAGE  TO  SOLVE  THE  PTCN0474 

SYSTEM.  THE  USER  MAY  WISH  TO  REPLACE  THIS  ROUTINE  WITH  ONE  MORE  PTCN047S 

SUITED  TO  HIS  PROBI.EM.  PTCN047A 

PTCN0477 

PTCN0478 

USER  SUPPLIED  SUBROUTINES  PTCN0479 

PTCN0480 

PTCN0481 

FCTNINVAR.X.FX)  PTCN0482 

PTCN04B.3 

EVAI.UATES  THE  NVAR-1  COMPONENT  FUNCTION  FX  GIVEN  X  AN  NVAR  COMPONENT  PTCN0484 

VECTOR.  THIS  FUNCTION  DESCRIBES  THE  NONLINEAR  SYSTEM.  THE  AUGMENTINGPTCM048S 

EQUATION  IS  HANDLED  BY  THE  CONTINUATION  PACKAGE.  PTCN0486 

PTCN0487 

FPR IME( NVAR.X. FPRYM)  PTCN048R 

PTCN0489 

EVALUATES  THE  NVAR-.1  BY  NVAR  JACOBIAN  MATRIX  FPRYM  (X)  PTCND490 
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AT  X  AND  STORFS  IT  IN  THF  NVAR  BY  NUAR  ARRAY  FRRYM» 

SO  THAT  FPRYHIIfJ)  CONTAINS  (BFCO  <!)/  BX  (J)). 

THE  LAST  ROM  OF  FPRYN  (FOR  THE  A(JGHENTIN(5  EQUATION)  IS  INSERTED 
BY  THE  ROUTINE  SOlUE. 


LINPAK  ROUTINES  USED 


LINPAK  reference: 

I.INPACK  USER'S  GUIDEj 

J  J  DONGARRAt  J  R  BUNCH •  C  B  HOLER  ANB  6  U  STEUART> 
SOCIETY  FOR  INDUSTRIAL  AND  APPLIED  HATHEMATICSf 
PHILADELPHIA!  1979. 


ISAHAXINfSXfINCX) 

INTEGER  FUNCTION  RETURNS  THE  POSITION  OF  LARGEST  ELEMENT  OF  SX 

SAXPYINiSAiSXiINCXiSYtINCY) 

SETS  VECTOR  SY<I)  *  SA*SX<I)+SY(I) 

SCOPY  ( N  >  SX  f  INCX  •  SY !  INCY ) 

SETS  SY(I)=SX(I) 

SnOT<N»SX.INCX!SrfINCY) 

SnOT  *  SIIH<I=1  TO  N)  SX(I)*SY(I) 

SNRH2(N>SX>INCX) 

SNRH2  »  EUCLIDEAN  NORN  OF  SX(I) 

**NOTE**  SNRH2  HAS  MACHINE  DEPENDENT  CUTOFF  CONSTANTS 

SSCAl.(NiSAfSX!lNCX) 

SETS  SX<I)*SA*SX(I) 

SGEFA(A>LDAfNfIPVT»JNFO) 

FACTORS  MATRIX  A  WHOSE  LEADING  DIMENSION  WAS  DECLARED  AS  LDA 
AND  WHOSE  ACTUAL  USED  DIMENSION  IS  N»  SETS  UP  PIVOT  INFORMATION 
IN  VECTOR  IPVT  AND  WARNS  OF  ZERO  PIVOTS. 


SOESl.IAfLDAfNfIPVTfBfJOB) 

ACCEPTS  OUTPUT  OF  S6EFA»  AND  A  RIGHT  HAND  SIDE  B,  AMD  SOLVES 
SYSTEM  A«X«B>  RETURNING  X  IN  B.  FOR  MODIFIED  NEWTON'S  METHODr  ONCE 
MATRIX  IS  FACTORED  BY  SGEFAt  ONLY  SGESl  IS  CALLED  FOR  SUCCESSIVE 
RIGHT  HAND  SIDES 


LABELED  COMMON  BLOCKS 


/COUNTl/  COUNTS  NUMBER  OF  CAI.LS  FROM  ...  TO  ...  AS  FOLLOWS: 
ICRSl.  »  CORECT  TO  SOLVE 

ITNSL  »  TANGNT  TO  SOLVE 

NSTCR  «  PITCON  TO  CORECT  FOR  IMPROVED  STARTING  POINT 
NCNCR  •  PITCON  TO  CORECT  FOR  CONTINUATION  POINT 
NTRCR  »  PITCON  TO  CORECT  FOR  TARGET  POINTS 

NLMCR  »  PITCON  TO  CORECT  FOR  LIMIT  POINT 

NLMRT  ■  PITCON  TO  ROOT  FOR  LIMIT  POINT 


NOTE  THAT  NSTCR*  NCNCR*  NTRCR*  Nt.MCR  AND  NLMRT  COUNT  THE  NUMBER 
OF  ITERATIVE  STEPS  (NEWTON  OR  R(M)TFINDING)  AND  NOT  JUST  THE  NUMBER 
OF  SUBROUTINE  CALLS. 

/COUNTS/  KEEPS  PERFORMANCE  AND  WORK  STATISTICS 
IFEVAL  »  NUMBFR  OF  CALl.S  TO  FCTN 

IPEVAL  «  NUMBER  OF  CALl.S  TO  FPRIME 

ISOLVE  »  NUMBFR  OF  CALLS  TO  SOLVE 

NRED  »  NUMBER  OF  STEPSIZE  REDUCTIONS  MADE  BEFORE  PREDICTOR 


PTCN049t 

PTCN049? 

PTCN0493 

PTCN0494 

PTrN0495 

PTCN049A 

PTCN0497 

PTCN0498 

PTCN0499 

PTCN0500 

PTCNOSO: 

PTCN0S02 

PTCN0503 

PTCN0504 

PTCN0S05 

PTCNOSOA 

PTCN0S07 

PTCN0508 

PTCN0509 

PTCN0510 

PTCNOSll 

PTCN0512 

PTCN0513 

PTCN0514 

PTCN0515 

PTCN051A 

PTCN0517 

PTCN0S18 

PTCN0519 

PTCN0520 

PTCN0521 

PTCN0522 

PTCM0523 

PTCN0524 

PTCN052S 

PTCN0S26 

PTCN0527 

PTrj(0528 

PTCM0529 

PTCN0530 

PTCN0531 

PTCN0S32 

PTCN05.13 

PTCN0534 

PTCN0S35 

PTCN0536 

PTCN0S37 

PTCN0S38 

PTCN0539 

PTCM0S40 

PTCN0541 

PTCN0542 

PTCM0543 

PTCN0544 

PTCNOS45 

PTCN0S4& 

PTCN0547 

PTrJ«0548 

PTCN0549 

PTCNO5S0 

PTCN05S1 

PTCN0552 

PTCN05B3 

PTCN05S4 

PTCN0S55 

PTCNOSSA 

PTCN05S7 

PTr.NOf.S3 

PTCN0S.S9 

ptcnosao 
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POINT  CONVERGED  TO  THE  NEW  CONTTNIIATTUn  POINT.  PTr,NO',ftl 

NRDSUN  =  TOTAL  NUMBER  OP  STEPSIZE  REDUCTIONS  PTCNOfiZ 

KN  =  NUMBER  OF  CORRECTOR  ITERATION  STEPS  TAKEN  IN  rtOST  RECENT  PTCNOS6S' 
CALL  TO  CORECT.  PTCN0r.b4 

KNSUM  »  TOTAL  NUMBER  OF  CORRECTOR  ITERATION  STEPS.  PTCN05A5 

PTCNOSAo 


/OUTPUT/  PTCNOSA7 

IWRITE  =  USER  ACCESSIBLE  OUTPUT  INDICATOR.  PTCNO^^S 

lURITF-sO.  NO  OUTPUT  PRINTED  BY  PITCON.  PTr,N0S<S9 

IMRITEalf  ERROR  MESSAGES  PRINTED  BY  PITCON.  PTCNOS70 

IHRITEa?»  CERTAIN  OUTPUT  HILL  BE  PRINTED  BY  PITCON.  PTr,N0571 

PTCN0S72 

/POINT/  CONTAINS  DATA  ABOUT  THE  SOLUTION  CURVE  PTCM057S 

DETA  »  BINARY  MANTISSA  OF  THE  DETERMINANT  OF  THE  AUGMENTED  JACOBIANPTCNOS74 

lEXP  a  BINARY  EXPONENT  OF  THE  DETERMINANT  OF  THE  AUGMENTED  JAC0BIANPTCN0S7S 

CURVCF  =  ESTIMATED  CURVATURE  BETWEEN  XC  AND  XF.  PTCN0576 

CORDIS  a  NORM  OF  THE  CORRECTOR  STEP  FROM  PREDICTED  POINT  TO  CORRECTEDPTCNOS?? 

POINTf  USING  MAXIMUM  ABSOLUTE  VALUE  AS  THE  NORM.  PTCN0S78 

THIS  QUANTITY  IS  MODIFIED  TO  AN  'OPTIMAL'  VALUE.  PTCN0S79 

ALPHLC  a  angle  BETWEEN  OLD  AND  NEW  TANGENTS  Tl.  AND  TC  PTCN05G0 

HSECLC  a  euclidean  MORN  OF  SECANT  BETWEEN  XL  AND  XC.  PTCN0S81 

FNRMXF  a  maximum  NORM  OF  FUNCTION  VALUE  AT  MEW  CONTINUATION  POINT.  PTCN0582 

PTCN0S8S 


/TOL/ 

EPMACHa  SMALLEST  NUMBER  SUCH  THAT  1 .OtEPMACH.GT.EPMACH 

..S*BETA«<1-TAU)  FOR  ROUNDEDf  TAU-DIGIT  ARITHMETIC 
BASE  BETA.  TWICE  THIS  VALUE  FOR  TRUNCATED  ARITHMETIC. 
THIS  IS  THE  RELATIVE  MACHINE  PRECISION. 

EPMACHa2**<-77)  FOR  DEC-10. 

EPSATEa  8BEPMACH 

EPSQRTa  SQUARE  ROOT  OF  EPMACH 


PROGRAMHING  NOTES 


THE  USER  MUST  - 

1.  WRITE  SUBROUTINES 

SUPPl.Y  A  CAl.l.INn  PROGRAM#  AND  THE  TWO  ROUTINES  FCTN  AND  FPRIME 
AS  DESCRIBED  ABOVE. 

2.  SET  STORAGE  AREAS 

DECLARE  A  REAL  VECTOR  RWORK  OF  SIZE  ISIZE#  ISIZE.GE.NVARXNVARfS) 
AND  AN  INTEGER  VECTOR  IPVT  OF  SIZE  NVAR. 

S.  PASS  CERTAIN  NON-DEFAULTABI.E  VAI.UES 

PASS  NVAR  GREATER  THAN  ZERO#  ISIZE.GE.NVARtlNVARfS) 

AND  SET  IRETaO#  KSTEP»-1  OR  KSTEP»0  ON  FIRST  CALI. 

FOR  A  NEW  PROBIEM. 

THE  USER  SHOUI.D  - 

1.  STORE  A  STARTING  POINT  XR  IN  THE  FIRST  NVAR  LOCATIONS  OF  RWORK 
BEFORE  CAl.LING  PITCON. 


PTCN0584 

PTCN058S 

PTCNOSSo 

PTCN0587 

PTCN0588 

PTCNOf.89 

PTCN0590 

PTCN0591 

PTCNOS92 

PTCNQS93 

PTCN0594 

PTnN0595 

PTCN0S96 

PTCN0597 

PTnN0S93 

PTCN0599 

PTCNOAOO 

PTCMOAOl 

PTCN0A02 

PTCN0A03 

PTCH0604 

PTCN060S 

PTCN0A06 

PTCN0AO7 

PTCN0A08 

PTCN0A09 

PTCNOAIO 

PTCN06U 

PTCN0A.t2 

PTCN0A13 

PTCNOA14 

PTCNOAIS 


IF  SUCH  A  VALUE  IS  NOT  GIVEN#  THE  CODE  NAY  BE  UNABI.E  TO  PRODUCE  ONE. 

2,  CAREFULLY  MONITOR  THE  VAl.UE  OF  IRET  ST  ’  AT  ANY  SERIOUS  ERROR 
IS  CAUGHT  BEFORE  ANOTHER  CAI.L  IS  MABE  TO  N. 

3.  CHOOSE  A  VALUE  OF  IMOD  FOR  THE  TYPE  OF  CORRECTOR  PROCESS  TO 
BE  USED. 


THE  USER  MAY  - 

1.  MONITOR  THE  PASSING  OF  BIFURCATION  POINTS  BY  SAVING  THE  OLD 
VALUF  OF  DETA  AND  COMPARING  IT  TO  THE  CURRENT  VAt.UE.  IF  THERE 
IS  A  CHANGE  IN  SIGN.  THEN  A  BIEIWCATTON  POINT  HAS  BEFH  PASSED. 


PTCNOftte 
PTnNOA17 
PTCN0A18 
PTCMOAl 9 
PTCNOAZO 
PTCNOA21 
PTCNOA2r 
PTCN0A23 
PTr,NOA24 
PTCNOAZS 
PTCNOA2A 
PTCMOA27 
PTCMOA2S 
PTCNOA29 
PTCNOASO 
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C  2.  ACCESS  THE  COHHON  BLOCKS  /COIIMTl/  AND  /C0tJNT2/  TO  KEEP  TRACK 
C  OF  THE  AMOUNT  OF  WORK  DONE. 

C 

C  3.  MONITOR  THE  COMMON  Bl OCK  /POINT/  FOR  INFORMATION  ABOUT  THE 
C  SOLUTION  CURVE. 

C 

C  4,  AT  ANY  TIMEj  RESET  THE  CODE  BY  PASSINO  IN  RSTEP=-1  OR  KSTEP=0. 

C  THIS  ALLOUS  THE  USER  TO  CHANRE  STEPSIZE.  DIRECTION  OF  CONTINUATIONt 
C  ERROR  CONTROLS*  OR  OTHER  PARAMETERS.  IT  ALSO  ENABLES  THE  USER 
C  TO  RUN  UNRELATED  PROBl.EMS  OF  DIFFERENT  SIZES  OR  ERROR  CONTROLS 
C  DURING  A  SINGLE  PROGRAM  EXECUTION. 

C 

c 

C  THIS  SUBROUTINE  IS  CALl.ED  BY 
C  USER  MAIN  PROGRAM 
C  AND  CALLS 

C  CORECT 

C  ROOT 

C  TANGNT 

C  FORTRAN  ABS 

C  FORTRAN  ACOS 

C  FORTRAN  ALOG 

C  FORTRAN  AMAXI 

C  FORTRAN  AMINl 

C  FORTRAN  DBLE 


PTCN0A31 

PTCN0632 

PTCN0433 

PTCN0A34 

PTCN0635 

PTCN063<S 

PTCNOA37 

PTCNOA3S 

PTCN0639 

PTCN0A40 

PTCN0<S4J 

PTCN0A42 

PTCN0A43 

PTCN0A44 

PTCN0A45 

PTCNOA46 

PTCN0c47 

PTCN0A48 

PTCM0<.49 

PTCN06S0 

PTCN0651 

PTCN04S2 

PTCN0A53 

PTCN0A54 

PTCN0655 


C  FORTRAN  FLOAT  PTCN0AS6 
C  FORTRAN  SIGN  PTCNOA57 
C  FORTRAN  SIN  PTCNOASO 
C  FORTRAN  SNGL  PTCN0AS9 
C  FORTRAN  SORT  PTCM06A0 
C  LINPAK  ISAMAX  PTCNOAaI 
C  LINPAK  SAXPY  PTCNOAA2 
C  l.INPAK  SCOPY  PTCN0663 
C  LINPAK  SNRN2  PTCN0AA4 
C  LINPAK  SSCAl.  PTCN0A6S 
C  PTCN0666 


INTEGER  IPVKNVAR) 

REAL.  RUORK(ISIZE) 

REAI.  URGE<8)>AC0F<12) 

DOUBLE  PRECISION  DTLIPC*DTCIPC*DAD,iUS*COSALF 
COMMON  /COUNTl/  ICRSL.ITNSLfNSTCR»NCNCR»NTRCR*NLMCR*NLMRT 
COMMON  /C0UNT2/  IFEVAl.*IPEVAl.»ISOl.VE*NRED*NRDSUN*KN*KMSUM 
COMMON  /OUTPUT/  lURITE 

COMMON  /POINT/  BETA . lEXP tCURVCF .CORDIS* Al.PHLC . HSECLC . FNRMXF 
COMMON  /TOl/  EPMACH.EPSATE.EPSORT 
DATA  I DONE  /O/ 

DATA  TENMl  /O.l/ 

DATA  TENM2  /O.Ol/ 

DATA  TENN3  /O.OOl/ 

DATA  URGE  / 

1  ,87351 lSE+00.  .1531947E+00.  .31918I5E-0I .  ,333994AE-I0. 

2  .4A77788E+00.  .A970123E-03.  ,19808A3E-05.  . tl22789E-08/ 


DATA  ACOF  / 

1  ,9043128E+00.- 

2  .851A099E+00.- 

3  .10400AlEf01. 


.7075A75E+00.-.4AA7383Ef01 .-.3A77482Et01 . 
.  1953n9E+00 .  - .  4830A3AE+01 .  - .  9770528E+00 . 
.3793395E-01.  .1042177E+01 .  .445070AE-0.t/ 


PTCN0AA8 

PTCN0AA9 

PTCN0A70 

PTCN0A71 

PTCN0A72 

PTCN0A73 

PTCN0A74 

PTCN0A75 

PTCN0A76 

PTCN0A77 

PTCNOA78 

PTCN0A79 

PTCN0A80 

PTCN0A81 

PTCN0A82 

PTCN0A83 

PTCN0A84 

PTCN0A85 

PTCN0A86 

PTCN0A87 

PTCN0A88 

PTCN0A89 


PTCN0691 


1.  PREPARATIONS.  PTCN0A92 

ON  FIRST  CALL  FOR  THIS  PROBLEM.  INITIALIZE  COUNTERS  AND  VARIABLES.  PTCM0A93 
CHECK  USER  INFORMATION  AND  SET  DEFAUI.TS.  AND  IF  <KSTEP.EO,-J > .  PTCN0694 

CHECK  NORM  OF  F(XR)  AND  CORRECT  XR  IF  NECESSARY.  PTnN0A95 

ON  EACH  CALL.  IF  INPUT  IRET  HAS  NONFATAL  VALUE.  RESET  IRET  PTCN0A96 

SO  THAT  CONTINUATION  LOOP  PICKS  UP  WHERE  IT  WAS  HALTED.  PTCN0A97 

PTrN0A98 

cmmmm««x«««««x«««sxxxx«»t«««x«x«tx»««*x««x««t*««xf««iii««xttt«4PTCN0A99 

C  PTCNOZOO 
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TERR^O 

rr  ( ri«ET.E0.-i>  tret*? 

IE  <IRET.E(J.-2.0R. rRET.E0.-3.0R.IRET.E0.-4)  ]RET*1 
IF  ( tRET,Ea.-5.(JR.TRET.Ea.-<i>  IRET*0 


IF  CODE  UA!}  CAELEO  ACATH  AFTER  FATAI.  *IALUE  OF  IRET» 

THEN  RETURN  WITH  ERROR  VALUE  IRET='10< 

TF  (IRET.LT.O)  GO  TO  440 

PERFORM  ONE-TIHE  ONLY  INITIALIZATIONS 

IF  (II10NE.NE.O)  GO  TO  10 

SET  THE  MACHINE  DEPENDENT  VARIABLE  EPMACH*  THE  SMALLEST  NUMBER 
50  THAT  <t .0+EPMACH,CT.l ,0) 

FOR  DEC  POP- 10  IN  SINGLE  PRECISION: 

EPMACH’'7.450!5804E-P 

FOR  IBM  360  OR  370  IN  SHORT  (SINOLE)  PRECISION: 
EPHACH*9.53674E-7 

FOR  COC  6600  OR  7400  IN  SINGLE  PRECISION: 

EPMACH«7.105427406t-15 

SET  EPSATE*3«EPNACH>  EPSQRT*SORT(EPMACH) 

EPSATE’S.OBEPMACH 

EPSORT*SORT(EPMACH) 

ALFMIN>2.0«AC0S(1 .O-EPMACH) 

IF<KSTEP,LT,-1 .OR.KSIEP.Gr.OlKSTEP*-! 

KSTEP0»-2 
I  DONE* 1 

PERFORM  INITIALIZATIONS  AND  CHECKS  FOR  NEW  PROBLEM  ONLY 

to  IF  <KST£P,OT.O)  GO  TO  30 

1F<KSTEP0.EQ.-1. AND. KSTEP.EO. 0)00  TO  30 

IF  (NMAR.LE.t)  GO  TO  440 

IF  (ISIZE.LT.<NVAR)«<NVAR45))  GO  TO  440 

IXR»1 

JXR*0 

IXC*IXR+NVAR 

JXCaJXR+NVAR 

TXFalXCfNVAR 

JXF»JXC4NVAR 

TTL*IXF+NVAR 

JTLaJXFFNVAR 

TTC*irLFNMAR 

JTCaJTLFNVAR 

IFPatTC+NMAR 

JFP»JTC+NVAR 

nETA*0.0 

TC1PC*0.0 

CORDlSaO.O 

CURVCFaO.O 

HSECLLaO.O 

HSECLC«0.0 

XTIO«0.0 

ITO-O 

NEON-NVAR-1 

NRED-O 

KNSIJN«0 

NR0SUN*0 

ICRSL*0 


PTCN0701 

PTCN0702 

PTCM0703 

PTCM0704 

PTCN0705 

PTCN07O6 

PTCN0707 

PTCN070G 

PTCN070V 

PTCN07t0 

PTCM07n 

PTCN0717 

PTCN0713 

PTCN0714 

PTCN0715 

PTCN0716 

PTCN0717 

PTCN0718 

PTCN0719 

PTCN0720 

PTCN0721 

PTCN0772 

PTCN0723 

PTCN0724 

PTCN0725 

PTCN0726 

PTCN0727 

PTCN072G 

PTCN0729 

PTCN0730 

PTCN0731 

PTCN0732 

PTCN0733 

PTCN0734 

PTCN073S 

PTCM073A 

PTCN0737 

PTCN0738 

PTCN0739 

PTCN0740 

PTCN0741 

PTCN0742 

PTCN074T 

PTCN0744 

PTCN0745 

PTCM0746 

PTCN0747 

PTCM074R 

PTCN0749 

PTCN0750 

PTCN075t 

PTCN0752 

PTCN07S3 

PTCN0754 

PTCN07SS 

PTCM075A 

PTCN0757 

PTCM0758 

PTCN07S9 

PTCN0760 

PTCN076t 

PTCM0762 

PTCM0763 

PTr.N0764 

PTCM0765 

PTCN076A 

PTCN07A7 

PTCN0760 

PTCN0769 

PTCM0770 
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NSTCR^O  PTCN077J 

PTCN077? 

MTRCR«0  PTCM0773 

(♦l.m:R»0  PTCN0774 

MLNRTsO  PTCN077r. 

IF»:VAI.-0  PTCM077A 

IS0l.yE»0  PTCR0777 

TPFMAl.aO  PTCN077R 

IF  (HMAX.LE.0.0)  HMAX=SaRT < FLOAT ( NVAR > )  PTCM0779 

IF  (HHIM.LFKPSURT)  HHtR=£PSORT  PTCN0780 

HDFFa.5*(HMAX+HHIM)  PTCR0781 

IF  (HFACT.LF-t .0)  HFACT-3.0  PTCN0787 

HRED-1 .0/HFACT  PTCM0783 

IF  (ABSERR.LE.0.0)  ABSERR^EPSORT  PTi:N0784 

IF  (RFLERR.LF.0.0)  RELERR=EPSnRT  PTrM0785 

IF  (IPC.LE.O.IJR.  IPC.iiT^RMAU)  IPC=«VAR  PTCN078A 

IF  (LIH.LT.O.OR.LIH.GT.NVAR)  LIH»0  PTCN0787 

IF  (H.ER.O.O)  H>HDEF  PTCN0788 

DI»-St6M<l .O.H)  PTCH0789 

HbABSIH)  PTCN0790 

PTCH0791 

IF  (KSTEP.LT.O)  CHECK  NORH  OF  F<XR)  AT  STARTING  POINT*  PTCN0792 

IF  ACCEPTABLE*  RETURN  rHNEOIATELY  WITH  KSTEP^O*  PTCN0793 

OTHERWISE  APPLY  NEWTON'S  HETHODt  HOLDING  VALUE  OF  PTCN079A 

TPC-TH  CUNPONCNT  FIXED.  PTCH079S 

PTCN079A 

IF  (KSTEP.GE.O)  60  TO  20  PTCH0797 

CALL  CORECT(NUAR>RWORK ( IXR) * IPC* RWORKI ITL  >  t lERR* INOD*RHORK( IFP) *  PTCN0798 
1  lPVT»A8«ERR*RELERR»XSTEP.N{:0N»FNttHXF)  PTCN0799 

NSTCR«NSTCR4KN  PTCN0800 

PTCN0801 

IF  NO  ACCEPTABLE  POINT  FOUND*  ERROR  RETURN  PTCN0802 

PTCN0803 

IF  (lERR.NE.O)  GO  TO  400  PTCN0804 

KST£PO«-t  PTCN0805 

KSTEPaO  PTCN080A 

HTANCF-H  PTCN0807 

GO  TO  340  PTCN080R 

20  IF  (KSTEP.EO.O)  CALL  SCOPY<H«;aR,RHORK( IXR) *  1 /RWORKI IXC) *  1 )  PTCN0809 

IF  (KSTEP.EO.O)  CALL  SCOPY(NMAR*RUORK( IXR) *  1 *RHORK( IXF > *  1 )  PTCN0810 

C  PTCNOOll 

C  PTCN0813 

C  2.  TARGET  POINT  CHECK.  IF  (IT.NE.O)  TARGET  POINTS  ARE  SOUGHT.  PTCN0814 

C  CHECK  TO  SEE  IF  TARGET  CONPONENT  IT  HAS  VALUE  XIT  I.YINO  PTCH0815 

C  BETWEEN  XC(IT)  AND  XF<IT).  IF  SO*  GET  LINEARLY  INTERPOLATED  PTCN08U 

r  STARTING  POINT.  ANO  USE  NEWTON'S  NETHOD  TO  GET  TARGET  POINT  PTCN08.I7 

C  PTCN081S 

CXXXXXXXXXXXXXXXtXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXtXtXXXXXXXXXXtXXXXtXXXXXXPTCN0819 


2.  TARGET  POINT  CHECK.  IF  (IT.NE.O)  TARGET  POINTS  ARE  SOUGHT. 
CHECK  TO  SEE  IF  TARGET  CONPONENT  IT  HAS  VALUE  XIT  LYING 
BETWEEN  XC(IT)  ANO  XF(IT).  IF  SO*  GET  LINEARLY  INTERPOLATED 
STARTING  POINT*  ANO  USE  NEWTON'S  NETHOD  TO  GET  TARGET  POINT 


30  IF(ir.LT.0.0R.TT.GT.HVAR)IT«0 
IF  (IT.EP.O)  GO  TO  40 

IF  (IRET.EO.t.AND.XIT.EO.XITO.ANO.IT.EO.irO)  GO  TO  40 
XCIT>RHORK(JXCtIT) 

XFTTaRHORK(JXFHT) 

IF  ((XIT.LT.XCIT).AND.(XIT.IT.XFIT))  GO  10  40 
IF  ((XTT.GT.XCIT).ANO.(XCT.QT.XFIT))  GO  TO  40 
DEL»XFIT-XCIT 


PTCN0820 

PTCN0821 

PTCN0822 

PTCN0823 

PTCH0824 

PTCH0825 

PTCN082A 

PTCN0827 

PTCN0828 


RATsO.O 

IF  (ABS(DEL),GT,EP80RT)  RAT»<XIT-XCIT)/nEL 
CALL  SCOPY ( NVAR  * RWORK ( IXF ) *  t  * RWORK ( IXR ) *  1 ) 

CALL  SAXP Y ( NVAR  *  - 1 . 0  *  RWORK ( IXC ) *  1  *  RWORK ( I XR ) *  1 ) 

CALL  SSCAL ( NVAR  *  RAT  *  RWORK ( IXR ) / 1 ) 

CALL  SAXPY(NVAR*1 .0*RW0RK(IXC)*l*RU0RK(IXR>*t> 

RWORK (IXR LIT) aXIT 

CALL  CORECT(NVAR*RWORK( IXR) * IT*RWORK( ITL) * lERR* INOD.RWORK ( IFP) * 
t  [PVT  *  ABSERR  *  RELERR  *  XSTEP*  NEON.FNRN ) 

NTRCR-NTRCR+KN 

iro»iT 

XTTn«XIT 


PTCN0829 

PTCN0830 

PTCN0831 

PTCN0B32 

PTCN0833 

PTCN0834 

PTCN0835 

PTCN083A 

PTCN0837 

PTCN0838 

PTCN0S39 

PTCH0840 


non 
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TF  (IFRR.FQ.O)  GO  TO  720  PTCH0B41 

IF  <  IFRR.irO.-  t )  00  TO  770  PTCN0847 

IF  (rrRR.E(J.-2)  GO  TO  3AO  PTr,N0843 

IK  <  00  TO  780  PTCM0844 

C  PTCNOGAS 

C  PTCN0847 

C  S,  T8NCEHT  AND  I.OCAI.  CONTINOATTON  PARAMETER  CALCULATION.  UNLESS  THEPTCN0848 

C  TANGENT  AND  LIMIT  POINT  CALCULATIONS  MERE  Al.REAOY  PERFORMED  (BECAUSE  PTCM0849 
C  THE  LOOP  WAS  INTERRUPTED  FOR  A  LIMIT  POINT)*  SET  UP  AND  SOLVE  PTCN0850 

C  THE  EOlJATION  FOR  THE  TANGENT  VECTOR.  FORCE  THE  TANGENT  VECTOR  TO  PE  PTCN08SI 

C  OF  UNIT  LENGTH*  AND  KORC.E  THE  IPL -COMPONENT  TO  HAVE  THE  SAME  SIGN  AS  PTCN08S2 
C  THE  IPl-TH  COMPONENT  OF  THE  PREVIOUS  TANGENT  VECTOR.  OR  (ON  FIRST  PTCN08SS 
C  STEP)  THE  SAME  SIGN  AS  THE  USER  INPUT  DIRECTION  DIR.  SET  THE  LOCAL  PTCN0854 
C  PARAMETER  IPC  TO  THE  1.0CATI0H  OF  THE  LARGEST  COMPOHEMT  PTCM08S5 

C  OF  THE  TANGENT  VECTOR.  UNLESS  A  LTNTT  POINT  IN  THAT  DIRECTION  PTCN08S6 

C  APPEARS  TO  PE  APPROACHING  AND  ANOTHER  CHOICE  IS  AVAIL ARLE.  PTCM0857 

C  PTCN08S8 

CXXXXtXXXXXktttXXttXttXXtXXXXXtXtXXtXXXtXtXtXXtXtXXXXXXXtXtXtXXtXXttttXXPTCHOSS? 


40  IF  <IRET.HF.2)  00  TO  50 
IRET=0 
GO  TO  UO 

STORE  OLD  TANGENT  TH  TL.  COMPUTE  NEW  TANGENT  EOR  XC 


PTCN08A0 

PTCN08<>1 

PTCN08A2 

PTCN08A.7 

PTCN08A4 

PTCN08A5 

PTCN08A6 

PTCN08A7 

PTCN08A8 

PTCN08A9 


50  IPl»IPC  PTCN08A7 

IF ( KSTEP . 6T . 0 ) CALL  SCOP Y ( NVAR » RMORK ( I TC ) . 1 . RHORK ( I TL ) / 1 )  PTCN08A8 

ICAI.L*!  PTCN08A9 

CALL  TANGNT( NVAR. RHORK ( fXF) . tPC.RMORKI ITC) > IRET. ICALL.RHORKI IFP).  PTCN0870 
1  TPVT.NEON.OETA.IEXP)  PTCN0871 

IF  (TRET.E0.-2)  GO  TO  430  PTCN0872 

IF  <IRET.E0,-1)  GO  TO  410  PTCN0873 

PTCN0874 

SUPRQUTINF  TANGNT  RETURNED  IPC.  THE  LOCATION  OF  THE  LARGEST  r0HP0NEHTPTCN0875 


SUPRQUTINE  TANGNT  RETURNED  IPC.  THE  LOCATION  OF  THE  LARGEST  COHPI 
OF  THE  TANGENT  VECTOR.  THIS  MILL  8E  USS-D  FOR  THE  LOCAL 
PARAMETERI7ATI0N  OF  THE  CURVE  UNLESS  A  LIMIT  POINT  IN  IPC  SEEMS 
TO  BE  COMING.  TO  CHECK  THIS.  ME  COMPARE  TClPCt«rC( IPC)  AND  THF 
SECOND  LARGEST  COMPONENT  TC.iPC:»TC( JPC) .  IF  TCJPC  IS  NO  LESS 
THAN  .1  OF  TCTPC.  AND  ICIIPC)  IS  LARGER  THAN  TL(JPC). 

WHEREAS  TCdPC)  IS  LESS  THAN  TLIIPC).  HE  HILL  RESET  THE 
LOCAL  PARANETER  IPCI--JPC. 

TLIPL*TCTPC 
TCIPC--RHORK<  JTC+IPC) 

IPC- IPC 

IF(ABS(TCIPC).GE.ABS<RHORK(.JTL+IPC)))GO  TO  AO 

TF<TLIPL,EO.O.O)GOTO  AO 

RUORK( JTC+IPC)®0.0 

.JPC=-- 1 S AMAX  <  NVAR .  RMORK  (  I  TC )  >  1 ) 

TC  JPC=RHORK  ( JTC+.IPC ) 

RMORK <jrCKTPC)*TCTPC 

IF  (ABS<TCJPC).LT.TENM1XABS<TCIPC))  GO  TO  AO 
IF  (ABS<rC.JPC).LT,ABS<RUORK<mKJPC)))GOTO  AO 
IPCsJPC 

TF  (IMRITE.GE.2)MRtTE(A.AtO) 

60  TCIPl»RHORK<JTCtIPL) 

DTI.  TPCsBBLE  ( RHORK  ( JTL  t  IPC ) ) 

DETA*DETA/TCIPL 

ADJUST  SIGN  OF  TANGENT 

COMPARE  THE  SION  OF  TC(TPL)  MTTH  SIGN  OF  TL(TPL) 

(BUT  ON  FIRST  STEP.  COMPARE  SIGN  OF  TCdPL)  HITH  USER  INPUT  DIR) 
IF  THESE  SIGNS  DIFFER.  CHANGE  THE  SION  OF  TC  TO  FORCE  AGREEMENT. 
AND  THE  SIGN  OF  DETA. 

THEN  RECORD  DIRJ*  SIGN  OF  OETERNINANr  »  SIGN(DETA). 

STLIPI.*0,IR 

IF  (TLIPL.ME.0.0)  STLIPL=SIGM( 1 .0. TLIPl ) 

IF  (SIGNd  .O.TCfPL)  .EO.STLIPL)  GO  TO  70 


PTCM087A 

PTCN0877 

PTCM0878 

PTCN0879 

PTCN0R80 

PTCN08B1 

PTCN0B82 

PTCN0883 

PTCN0884 

PTCN0885 

PTCN08BA 

PTCN08B7 

PrCN0888 

PTCM0889 

PTCN0890 

PTCN0891 

PTCN0897 

PTCN0893 

PTCN0894 

PTCN0895 

PTCN089A 

PTCN0897 

PTCN0898 

PTCN0899 

PTCH0900 

PTCN0901 

PTCN0902 

PTC.M0903 

PTCN0904 

PTCH0905 

PTrN09OA 

PTCM0907 

PrCN0908 

PTCN0909 

PTCN0910 
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CAI.L  SSCAKNVAR.-I  .Otf(UORK<  ITO  f }  > 

1>ETA*-DETA 

TniPl»-TCIPl. 

70  T(;TPC«RHl)RK<.JTt:+rPC) 

TC  JPC--RW0RK  <  JTC + JPC ) 

0TrTPC«08LEai:tPC) 

niRsSr6M<1.0fPPTA) 

IP  <i.TH.£U>0)  Hi)  10  80 
Tl.  LIMsTCLIH 
Ta.rH'«RUORK<.)TCRrH) 

COMPUTE  ALPHI.C.  THE  AMOLE  8ETMEEM  TAMGEMT  AT  XL  AMO  TAM6ENT  AT  XC 

AMP  HSECLC>  THE  EUCLIDEAN  NORM  OF  SECANT  FROM  XL  TO  XC. 

80  rF(KSTfP.LE.O)  GO  TO  140 
C08ALK-0.000 
DO  90  I«lfNVAR 

COSALF^COSALFf OBLE < RMORX ( JTCf [ ) ) BDBLE ( RUORK (  J  TL  M ) ) 

90  RUORK ( JXRP I ) «RUORK  <  JXFP I ) -RUORK ( JXCP I ) 

HSECLL=HSECLC 

HSECLC>SNRN2(NUARfRH0RK(IXR)>l > 

ALPHI.C^SNGL<CQSALF) 

IFIALPHLC.GT. 1 .0)ALPHIC«] .0 
IFIALPHLC.I-T.  -1 ,0)ALPHLC«-1 .0 
ALPHLC«ACOS<ALPHLC) 
rF<T«RrTE.CE.2»«PITE<4.550>AI.PHLC 


PTCN0911 
PTCN0912 
PTCN0913 
PTCN091A 
PTCM0915 
PTCHOVt A 
PTCN0917 
PTCM091R 
PTCN091 9 
PTCM0920 
PTCN0921 
PTCH0972 
PTCN0923 
PTCN0924 
PTCN0925 
PTCM0976 
PTCN0927 
PTCM0928 
PTCM0929 
PTCN09.T0 
PTCN0931 
PTCM0932 
PTCN0933 
PTCN0934 
PTCM0935 
PTCN0934 
PTCM0937 


4,  LIMIT  POINT  CHECK.  IF  (LIM.ME.O)  CHECK  TO  SEE  IF 
OLD  AND  NEU  TANGENTS  DIFFER  IN  SION  OF  LIM>TH  COMPONENT. 
IF  S0»  ATTEMPT  TO  COMPUTE  A  POINT  X»  8ETHEEM  XC  AMO  XF 
FOR  WHICH  TANGENT  COMPONENT  VANISHES 


IF  (LIM.LI.O.OR.KSTEP.IE.O)  GO  TO  160 
CHECK  FOR  LIMIT  INTERVAL 

IF  <SIGN(1.0«TCi.IM).Ea.STGN<1.0tTLLIM>)  GO  TO  160 
TEST  FOR  ACCEPTABLE  ENDPOINTS 
ATl.LN>ABS(TLLIN) 

IF  <ATLLM.GT,0,5»A8S£»R)  00  TO  110 

IF  XC  IS  LIMIT  POIMf*  TL  ALREADY  CONTAINS  TANGENT  AT  XC 

100  CALL  SCOPY(NVAR>RWORK<  (XOft/RHORXdXR)*!) 

60  TO  310 

110  ATCLIN»ABS<TCLIM) 

IF  (ATCLIH.GT.O.SBABSERR)  GO  TO  130 
120  CALL  SCOPY(NVAR»RUORK(TXF)»t»RHORX<rXR>/l) 

CALL  SCOPY  <  NVAR • RUORK  < 1 TC ) 1 1 f  RUORK  <  T  TL  >  f 1 ) 

GO  TO  310 

TEST  FOR  SMALL  INTERVAL 

130  XCL(H-RUORK^  :XCH  (H) 

XFL1N«RU0RK(JXFHIN) 

BCL«ABS<XFLTN-XCLIN) 

XAB8« AMAX 1 <  ABS ( XCL I N ) f  ABS  <  XFL I M ) > 

IF  <DEL.0r>(ABSERR4RELERR«XABS))  GO  TO  140 
IF  TATLLN.OT.ATCLIN)  00  TO  120 
GO  TO  100 

BEGIN  RUOT-FINOING  ITERATION  WITH  INTERVAL  (Od)  AND 
FUMCTIOU  VALUES  Tl.ILIM).  TC(LIM). 


PTCN0949 

PTCN09S0 


PTCN0951 

PTCN09S2 

PTCM0953 

PTCM09S4 

PTCM0955 

PTCN09S6 


PTCN0957 

PTCN09S8 

PTCM0959 

PTCN0960 

PTCM0961 

PTCN09A2 

PTCN0963 

PTCN0964 


PTCM0965 

PTCM09A6 


PTCM0967 

PTCN0968 


PTCM0969 

PTCM0970 

PTCN0971 

PTCM0972 

PTCN0973 

PTCN0974 

PTCM0975 

PTCN0976 


PTCN0977 


PTCN097R 

PTCM0979 

PTCN0980 
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140  K0«IMT*0 
AnO.O 
FA«TLLIM 
TSM^^raiH 
C«1 .0 
FCsTCl.IM 


PTCM0981 

PTCH0987 

PTCM0983 

PTCN0984 

PTCM0985 

PTi:;>J098A 


SkT  rp».fN  TO  THK  xhm:(  of  haxthoh  or»y  of  hfcaht 
(FXCFPT  THAT  TPLIH  HliST  NOT  EQUAL  t.IH) 

AND  8AUE  THE  SIGN  OF  THE  HAXTHUH  COHPOHEHT  TH  DIKLIH 
so  THAT  NEH  TANGENTS  HAY  RE  PROPERLY  SIGHED. 

CALL  SCOPY(NUAR»RUORK(IXF)fl rRUORK(JXR)>l ) 

CALI.  SAXPY  ( NUAR  i.  - 1  , 0 1 RMURX  (  TXC » ♦  1  »  RUORK  <  t  XR  > .  1  I 

RUORK(  JXRHIH)’°0.0 

TPI.  tH*ISAHAX(MUARtRUOWK(  TXX)  >  1 ) 

DIRLIHaSTGN(1.0tRUORK< JXR4IPL;H) > 

CALL  ROOTFINDER  FOR  APPROXIHATE  ROOT  SNt  SET  X^^SNtXF-K  L-SN)«XC 
CALL  CORRECTOR  TO  RETURN  TQ  CURVE  OH  LINE  X( TPI.  [H)=CONSTANTf 
COHPUTE  TANGENT  THERE t  AND  SET  FUNCTION  VALUE  TO  TANGENT (LIH) 

150  CALL  R00T(A»FA»SNtTSN»CfFC>K0UNT>IFLA6> 

NLHRT«NL HRTLl 
IF(  fFLAG.LI.-DGO  TO  350 
IFdFLAG.EP.-l  .OR.IFLAG.EO.OIGO  TO  310 
CALL  SCOP Y ( NUAR » RWORK ( IXF ) < 1 >  RUORK ( IXR  >  d ) 

CALL  SSCAL INVAR »SNfRHORK(]XR)>l) 

SCALER*1.0-SN 

CALL  S AXP Y < NVAR • SC ALER f RUORK < I XC ) > 1 > RUORK ( I XR ) > 1 ) 


PTCN0987 

PTCN09SB 

PTCN0989 

PTCN0990 

PTCN0991 

PTCN0992 

PTCN0993 

PTCN0994 

PTCN0995 

PTrN099A 

PTCN0997 

PTCN0998 

PTCN0999 

PTCNIOOO 

PTCNIOOI 

PTCN1002 

PTCN1003 

PTCN1004 

PTCNIOOS 

PTCHIOOA 

PTCN1007 

PTCN1008 

PTCN1009 

PTCNIOIO 


CALL  CORECT ( NVAR • RUORK < IXR ) > TPL IN ^ RUORK < I TL > • TERR • INOO* RUORK ITFP >  PTCN 1 A U 


PTCH1012 
PTCNlOtX 
PTCN1014 
PTCNtOlS 
PTCN10J4 
PTCN1017 
PTCN1018 
PTCM1019 
PTCN1020 
PTCN10?t 
PTCN1022 
PTCN1023 
PTCN1024 
PTCNt025 
PTCN102A 
PTCHt027 
PTCN1028 
PTCN1029 
PTCN1030 
PTCNlOXt 
PTCNt032 
PTCNt033 
PTCN1034 
PTCN 1035 

PTCNt037 

5.  STEP  LENGTH  COHPUTATION,  COHPUTE  8TEPLENGTH  HTAHCF  PTCN 1038 

PTCN1039 
PTCN1040 
PTCN1041 

LET  PTCN1042 

PTCH1043 

ALPHLC  «  HAXINUH  OF  ARCCOSITLf TO  AND  At.FHIN  '  2«ARCC0S(1-EPNACH)  PTCN1044 

HSECLC  -  NORN(XL-XC)  PTCN1045 

HSECLL  »  NORH<XL-XLL)  PTCN104* 

ABSIN  “  ABSIGINI ,5*ALPHLC>)  PTCM1047 

CURVLC  «  LAST  VALUE  OF  CURVCF  PTCN1048 

CURVCF  «  2XABSTN/HSECLC  PTCN1049 

CORDIS  *  (IPTIHI7ED  CORRECTOR  DISTANCE  TO  CURRENT  CONTINUATION  POINT.  PTCN1050 


1  . TPVT  >  ABSERR  >  REL  ERR • XSTPLN • NEON » F  NRH) 

NLNCR*NLNCR+KN 

IF  (lERR.NE.O)  GO  TO  350 

TCALL«1 

IPT»IPLIH 

CALL  rAN6Nr(NVARtRU0RK<rXR)«IPr>RU0RK(TTI.)>  IREIdCALL  . 

1  RU(IRK<IFP)»IPVT*NEQN>DETLIH>IEXLIH} 

IFdRET.LT.OIGO  TO  350 

AQ.JUST  THE  SIGN  OF  THE  TANGENT  SO  THAT  THE  TPLIN-TH  CONPONENT 
HAS  THE  SANE  SIGN  AS  THE  IPLIN-TH  CONPONENT  OF  THE  SECANT 

IF  (UIRLIN.NE.SIGNd.OfRUORKIJTLFIPlIH))) 

1  CALL  SSCAL  < NVAR »-1.0>  RUORK  dTDd) 

SEE  IF  UE  CAN  ACCEPT  THE  NEW  POINT  BECAUSE  TANGENKLIN)  IS  SHALL 
OR  NUST  'ACCEPT'  THE  POINT  BECAUSE  THE  VALUES  ARE  HOT  DECREASING 
RAPIDLY  ENOUGH*  OR  IF  UE  CAN  GO  ON. 

TSNOLDaTSN 

TSN>RUaRK(JTL-M.lN) 

IF  (ABSITSN) .LE.O.SXABSERR)  GO  TO  310 
GO  TO  150 


THE  FORHULAS  UNDERLYING  THE  AL60RITHH  ARE 
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BUT  CORDIS  FORCFD  70  LIE  BETWEEN  .ItHSEri-C  AND  HSFCl.C. 
UNLESS  (COROTS. ED. 0.0)7  BECAUSE  THE  BREOtCTEO  POINT  WAS 
INNEDIATELT  TTrCEPTED.  JN  SUCH  A  CASE^  SET  HTAHCE=HSECI.  C 
INSTEAD  OE  USTHU  FIRST  ESTIMATE  EOK  HTANCF. 


THEN 


CURUXE  -  CURUCK  +  HSECl.CtlCURUCF-CURVI.O/IHSECl.CfHSECLl  ) 

BUT  CIIRUXE  MUST  BE  GREATER  THAN  .OOlr  AND  A  SIMPLER  FORMULA  IS  USED 
IE  WE  DO  NOT  HAUE  DATA  AT  TWO  PRErUTOUS  POINTS. 

FIRST  ESTIMATE?  (UNLESS  (CORDIS, EQ. 0,0)  ) 


HTANCF 


SORT  ( i’BCOROIS/CURUXF  ) 


PTCM1051 
PTCN1052 
PTCNJOt^S 
PTCNIOSA 
PTCM1055 
PTCNtOSA 
PTCN1057 
PTCNIOSR 
PTCM1059 
PTCNI040 
PTCN)06) 
PTCN10A2 
PTCN106.T 
PTCN10A4 
PTCN10(i5 
PTCNIOAA 
PTCM10A7 
PTCN1068 
PTCN10A9 
PTCM1070 
PTCM1071 
PTr.Nl072 
PTCN1073 
PTCH1074 
PTCM1075 
PTCM107<4 
PTCN1077 
PTCN1078 

lEttttktttttlttBIlBBttBtXtttXttBttttBXBBBttBBtBtBBttBBBtBBttttBtBtttttXttPTCNlOT? 

PTCM1080 

CHECK  IF  DEFAULT  STEP  MUST  BE  USED! 

ON  FIRST  STFP*  HTANCF^H* 

IF  PREUIOUS  step' WAS  OF  SIZE  ZEROt  USE  STEPSTZE  H0EF»(HMTN+HMAX)/2 


ADJUSTED  UALUE: 

HTANCF  a  HTANCFtd.O  ♦  HTANCF»(  TC(  IPC  )-TL(  IPC ) )/(  ?»HSECI.  CkTCI  IPC  > ) ) 
READJUSTMENT  AND  TRUNCAT TONS.' 

IF  STEPSTZE  REDUCTION  OCCURRED  DURING  LAST  CORRECTOR  PROCESS. 

HTANCF  IS  FORCED  TO  BE  LESS  THAN  (HFACT-t ) »HSECLC/2 , 

HTANCF  MUST  LIE  DETWEEN  (HSECLC/HFACT)  AND  (HSECLCkHFACT ) , 

HTANCF  IS  ALWAYS  FORCED  TO  LIE  BETWEEN  HHIN  AND  HHAX. 


tAO  If  (KSrEP,8T,0.AND.HSECLC.GT,0.0)  00  TO  170 
HTANCF»HDEF 
IF(KSTEP.LE,0)HTAMCF»H 
GO  TO  190 

170  IF(ALPHLC,LT,ALFMIN)ALPHl.CaALFMIN 
APSTNaABS(SIN< .5«ALPHLC) ) 

CONPUTE  NEW  CURVATURE  DATA 

CURVLC»CURVCF 
CURWCFa2 . OBABS IN/H8ECLC 
CURVXF-CURVCF 
IF(HSECLL.NE.O.O) 

1  CURVXF«CURVCF+HSECLC*<CURVCE-CURVLC)/(HSECl.r.tHSECl.l. ) 
CURVXF«AMAX 1 ( CURVXF . TE NM3 ) 

IF(IWRITF.0E.?)MRITE<4.5A0) CURVCE . CURVXF 

IF  THE  CONVFRRENCE  DISTANCE  IS  ZEROt  SET  FIRST  ESTIMATE  TO  HSECLC. 
OTHERWISE t  TRUNCATE  CORDIS  TO  LIE  BETWEEN  .OttHSECLC  AND  NSECLC. 

HTANCF«HSECLC 

lE(CORDIS.EO.O.O)  GO  TO  180 
TEMP-=TENM2*HSECLC 
COROISaAMAXK CORDIS. TEMP) 

CORDISaAMINl (COROIStHSECLC) 

SET  HTAMCE.  THEN  ADJUST  FOR  CURVATURE  IN  CONTINUATION 
PARAMETER  DTRECTIONi  THEN  TRUNCATE 

HTANCF«SORT(2<0«C0R0rS/CURVXF> 

180  IF(NREO,6T.O)HrANCEaANINt(HrANCFt(HEACT~l,0)»HSECI.C»,5) 

DADJUSal.ODO+(l,ODO-DTLIPC/nTCTPC)»nm.E(  .5*HTANCF)/nBLE(HSEr.LC) 
HTANCF«NTANCF«SNGL ( DAD  JUS ) 

TEMPaHSECLCXHRED 
HTANCr-AMA?M  (HTANCF. lEMP) 

TFMPaMSECLCXHFACT 


PTCN1081 
PTCN1082 
PTCM1083 
PTCN1084 
PTCN1085 
PTCN1084 
PTCN1087 
PTCN1088 
PTCN1089 
PTCNJ090 
PTCN1091 
PTCNJ0V2 
PTCM1093 
PTCM1094 
PTCN1095 
PTCH1094 
PTCH1097 
PTCM1098 
PTCN1099 
PTCNllOO 
PTCHllOl 
PTCN1102 
PTCN1103 
PTCN1104 
PTCN110S 
PTCNllOA 
PTCNltOr 
PTCN1108 
PTCNlt09 
PTCNinO 
PTCHl 1 1 1 
PTCNllJZ 
PTCN1113 
PTCN1114 
PTCHlllS 
PTCMl  lU 
PTCN1117 
PTCMlllS 
PTCNI ( (9 
PTCNllZO 
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PTCNn21 
PTCNII?? 
PTCNn23 

C  PTCHtt74 


HTAMCF»AHIMi<HTAMrF.TEHP) 
HT AWCF-AMAX  t  < HTANCF , HM IN i 
HTAMCF^AHTMl ( HTAMCF i HMAX ) 


C 

C 

C 

C 

C 

C 

C 


A.  PREDICTION  AND  CORRECTION  STEPS.  USIMO  XRrXC+HTANCFXTC 
AS  STARTING  POINT?  CORRECT  XR  WITH  A  FOLE  OR  NOttTFIED  NEUTON 
ITERATION.  IF  CORECT  FAIl.Sf  REDUCE  STEPSI7E  USED  FOR  PREDICTOR 


POINT?  AND  TRY  AGAIN. 
FALLS  BELOW  HNIN. 


CORRECTION  WILL  ONLY  BE  ABANDONED  TF  STEPS  I ZE 


PrCN1126 

PTCNIIJ? 

PTCNn28 

PTCNlll’? 

PTCN1130 

PTCN1131 

PTCN11.T2 


CXXXXXXXXXXXXXXXXXXXXXkXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXtXXXXXXPTCN1133 


190 


200 

210 


PTCN1134 
PTCN113S 
PTCN1136 
PTCN1137 
PTCN1138 
PTCNl 139 
PTCN1140 
PTCNl 141 
PTCN1142 
PTCN1143 
PTCNl 144 
PTCNl 145 
PTCNl 144 
PTCNl 147 
PTCNl 148 
PTCN1149 
PTCNl 150 
PTCN1151 
PTCN1152 
PTCN1153 
PTCNl 154 
PTCNl 155 
PTCNtl54 
PTCNl 157 
PTCNl 158 

C  PTCN1159 

CXXXXXXXXXXtXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX>XXXXXXXXXXXXXXXXXXXXXPTCN1140 


KSTEPO»KSTEP 
KSTJ-P-KSTEPFl 
NREDaO 

CALL  SCOPY(NyARjRWORK( tXF > ? 1 .«WORK( IXR) ? 1 ) 

CALL  SAXP Y  ( NOAR  i NT ANCF ,  RUORK  <  I TC  >  r  1  > RWORK  ( I XR  >  ?  1 ) 

T F  < I MR  r IE . GE . 2 ) HR  t  TE ( 4 » 570  > HTANCF 
IF(TURITE.0E.3)URTTE(4«580)<RW0RK(JXR4I)>I=1 tNOAR) 

CALL  CORECT(NyAR>RUORKaXR).IFC>RWORK(TTL)?rERK?  [NOO?RUORKaFP >  ? 
1  IPUT  •  ADSFRR  f  REI.ERR  •  XSTEP .  NEON » FNRHXF ) 

NCNCR -NCHCK4KN 
IF  (lERR.EP.O)  60  TO  230 
IF  (lERR.EO.-t)  GO  TO  420 

NO  CONyERCENCE?  TRY  A  SNALL.ER  STEPS T7E 

htancf-hredxhtancf 

IF  <HTANCF.LT, HNIN)  60  TO  390 
NRED=N»E0+1 
IF  <IERR. EO. -2)  GO  TO 
GO  TO  200 

CALL  SAXPYINMAR?-] .0>RH0RK(IXC>?lrRW0RK(IXR>f 1 > 

CALI.  SSCAL  <  NOAR » HRED » RUORK  (  UR  )  » 1 ) 

CAL  L  SAXP Y ( NO AR  f 1 . 0  ?  RUORK ( I XC ) > 1 >  RUORK ( I XR  >  > 1 > 

GO  TO  210 


220 


220 


7.  SUCCESSFUL  STEP,  STORE  INFORHATION  BEFORE  RET)JRN, 
UPDATE  OLD  AND  CURRENT  CONTINUATION  POINTS. 

CONPIJU  CORDIS?  THE  SIZE  OF  THE  CORRECTOR  STEP,  COMPOTE 
A  FACTOR  THETA  WHICH  RESCALES  CORDIS  TO  A  VALUE  WHICH  WOULD 
CORRESPOND  TO  A  DESTKABL.S  NUNBER  OF  CORRECTOR  STEPS 
<4  FOR  FULL  NEWTON?  10  FOR  NODIFIED  NEWTON). 

SE.E  RE.^EReNCE  DEN  HEIJER  AND  RHETNBOLOT?  LOC,  CIT. 


PTnN1141 
PTCN1142 
PTCNl 143 
PTCNl 144 
PTCNl 145 
PTCNl 144 
PTCNl 147 
PTCNl  li48 
PTCNH49 


CXXXXXXXXXXXXXXXXXXXXXtXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXtXXXXXXXXPTCNtl70 


c 

c 

c 

c 


230  NR0SUN«NR»SUN4NRED 

IF(NRF0.NF.0.AND.IURITE.GE.2)URITE(4?590>NRED 


COMPUTE  CORRECTOR  STEP  * 
SET  CORDIS  ^  MAX  NORM  OF 


XC+HTANCrXTC-XF 
CORRECTOR  STEP 


C 

C 

C 

C 


CALL  SCOPY (NMAR? RWORK <  TXF) ? 1 /RUORK ( IXC) ?1 ) 

C AL  L  S AXPY  <  NV AR I - 1 . 0  ?  RUORK ( I XR ) ? 1 ?  RWORK  <  T  XC ) ?  I ) 
CALI.  SAXP  Y  <  NVAR  ?  HTANCF  *  RWORK  <  t  TC )  ?  t  ?  RUORK  ( t  XC  >  ?  1 ) 

I NAX- T  SANAX ( NVAR  ?  RUORK ( I XC ) ? 1 ) 

CORD TS »ABS ( RUORK  <  iXC 4 INAX ) ) 

IF<KN.FO.O)  CORDIS-O.O 

NODIFY  CORDIS  TO  A  VALUE  THAT  WOULD  CORRESPOND  TO  THE 
DESIRED  NUMBER  OF  CORRECTOR  STEPS 

rF^COR0TS.E«,O,0>GO  TO  300 

OMFGAsXSTFP/CORDTS 

THETA^O.O 


PTCN1171 
PTCMfl72 
PTCM1173 
PTCNl 174 
PTCNU75 
PTCMlt7A 
PTCNM77 
PTCNl 178 
PTCNl! 79 
PTCNH80 
PTCNn81 
PTCNl 182 
PTCNl 183 
PTCNl 184 
PTCNl 185 
PTCNl 184 
PTCNl 187 
PTCNl 188 
PTCN1189 
PTCNl 190 


oooonoooooo  ooo  oooo  ooo  nnn 
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IFdHOn.FQ.  1  )  GOTO  260 

FULL  NCUTON  METHOD  FOR  rORRECTOR  STEPS 

IF(KN,LE. J )THETA*8.0 
IF(KH.EH.4)THE rA=t»0 
IF<THETA. ME. 0.0)60  TO  290 
rF(KM.fjT.4)60  TO  240 
LK=4»KN-7 
THETA- t .0 

IF<0ME0A.GE.URGE(LK) )G0T0  290 
LST-l.K 

IF ( OMEGA. CE. URGE (LK4] ) )G0T0  250 
I.ST--LK+2 

IF^0NE6A.GE.URGE<LKi2) )G0T0  250 
THETA=8.0 
GOTO  290 
240  THETA=0.t25 

IF(KM.6E.7)  GOTO  290 
I.K-4«KN-U 

IF(OMEGA.LE.URGE(l.K)  )G0TQ  290 
LST-2*KM-.t 

250  THETA«AC0F<l.ST)4AC0F(LST4] )«AL06< OMEGA) 

GOTO  290 

MOOtFIEO  MEUTOM  METHOD  FOR  CORRECTOR  STEPS 

240  tF(KM,LE,l) THETA ~«.0 
IF<KM.EQ.10)THETA=1 .0 
IF(rHETA, ME. 0,0)00  TO  290 
EXPQN>>F1.0AT(KN-]  )/FLnAT(KN-10) 

AVOID  0VERF1.0U  OR  IIMDERFIOU  BY  ANTICIPATIMO 
CUTOFF  VAEUES  OF  THETA 

IF  (KN.GT.IO)  GO  TO  270 
IF  (8.0««EXP0N.GT.0MEGA>  THETA»8.0 
IF  (,125**EXP0N,I.T. OMEGA)  THETA*. 125 
IF  (THETA. ME. 0.0)  60  TO  290 
GO  TO  280 

270  IF  (G.OmXPOH.LT.  OMEGA)  THFTA*8.0 
IF  (,125»*EXP0M.GT, OMEGA)  THETA*. 125 
IF  (THETA. ME.O.O)  GO  TO  290 
280  EXP0M*1,0/EXP0M 

THETA*QME6A««EXP0N 
THETA*AMAX1( THETA70.125) 

THETA*AMIM1 (THETA. 8.0) 

SET  THE  MODIFIED  VALUE  OF  CORDIS 

290  CORDlSaTHETA«COROJS 

TF(  rURTTE.GE.2)URTTE(4.400)0MEGA/T.HETA. CORDIS 
300  CALL  SCOPY(NVARfRUORK( JXF).t.RUORK(IXC).l) 
CALL  SCOPY(MVAK.RUORK( IXR) • 1 .RUORK( IXF) . 1 ) 

60  TO  340 


RETURNS.  SET  VALUE  OF  IRET.  IF  AN  ERROR  OCCURRED.  PRINT 
A  MESSAGE  AS  HELL. 


RETURN  LIMIT  POINT 

310  IRET*2 
RETURN 

C 

C  RETURN  UTTH  TARGET  POINT 


PTCM1191 
PTi:N1192 
PTCM1193 
PTCN1194 
PTCNl J95 
PTCN1194 
PTCMH97 
PTCNl 198 
PTCNl 199 
PTCN1200 
PTCNl 201 
PTCNl 202 
PTCN1203 
PTCNl 204 
PTCNl 205 
PTCN1204 
PTCN1207 
.'»TCN1208 
PTCN1209 
PTCNl 210 
PTCNl 21 1 
PTCNl 21 2 
PTCN1213 
PTCN1214 
PTCNl 21 5 
PTCH1214 
PTCN1217 
PTCN1218 
PTCN1219 
PTCNl 220 
PTCN1221 
PTCN1222 
PTCN1223 
PTCN1224 
PTCNl 225 
PTCNl 224 
PTCN1227 
PTCNl 228 
PTCNl 229 
PTCN1230 
PTCN1231 
PTCN1232 
PTCN1233 
PTCN1234 
PTCNl 235 
PT.':n1234 
PICN1237 
PTCNl 238 
PTCN1239 
PTCN1240 
PTCNl 241 
PTCM1242 
PTCN1243 
PTCN1244 
PTCNl 245 
PTCN1244 
;PTCN1247 
PTCNl 248 
PTCN1249 
PTCNl 250 
PTCNl 251 
PTCNl 252 
PTCNl 253 
PTCM1254 
PTCNl 255 
PTCN1254 
PTCNl 257 
PTCNl 258 
PTCMi::59 
PTrN1240 


UCJC  CJOU 
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7^?.»  {K>:T=t 
RFTIIRM 

RETURN  UITH  CONTINUATION  POINT 

330  CALL  SCOPY(NVAR>RWORK(IXF)iliRUORK<IXR)f J > 

340  IRFT=0 
HsHTANCF 
RETURN 

error  RETURNS 
350  rRFr=-! 

IFdURITF.GE.  I  )URITE(^f4S0) 

RF  TURN 
3^0  IRFT=-2 

IF( rWRTTF.GF. t) WRITE <4»4A0) 

RETURN 
370  rREr--3 

IFdURITF.GE.I  )URITE(A«470) 

RETURN 
380  IRET*-4 

TFdURTTE>GE.l)WRTTE(Af480) 

RETltRN 
390  IRET:=-5 

IFdWRITE.GE.  l)URITE(A*490)HTANCFfHNIN 
RETURN 
400  IRET=-A 

IF(  rWRITE.GEd)URTTE(A>SOO) 

RETURN 
410  rRET^-7 

IFdURITE.GE<  1  )URJTE(AtS10) 

RETURN 
470  IRF.Ts-8 

If-dWRITE.GE,l  >WRtTE(A>520) 

RETURN 
430  TRET--9 

IFdURITE.GE.  1  )URITE(AfS30) 

RETURN 
440  IRET=-10 

IF  dWR  ITE .  GE .  1 )  WR  r  TE  <  A  »  540  )  NUAR » ISIZE 

rftlirn 

450  format <2AH0I.TNTT  POINT  FINDER  FAILED) 

440  FORMAT (50H0C0RECT.  SEEKING  TARGET  POINT f  TOOK  TOO  MANY  STEPS) 
470  FORMAT <50H0Cl)RECT»  SEEKING  TARGET»  CALLED  SOLOE  WHICH  FAILED) 
480  FORMAT (45H0C0RECT»  SEEKING  TARGET*  FAILED  WITH  RAD  STEP) 

490  FORMAT! lOHOSTEPSIZE  *Ft2.7*t5H  LESS  THAN  HMINyF12>7) 

500  FORMAT <47H0N0RM  OF  F(X)  IS  TOO  LARGE  ON  INITIAL  CALL) 

510  FORMAT (33H0S0LUE  FAILED  IN  CALL  FROM  TANGNT) 

520  FORMAT (33H0S0LWE  FAILED  IN  CALL  FROM  CORECT) 

530  FORMAT! 23HOTANGENT  UECTOR  IS  ZERO) 

540  FORMAT !2AH0UNACCEPTABLE  INPUT  NUAR>*I10f7H  ISI7E=fIJ0) 

550  FORMAT !34H  AMOLE  DETWEEN  OLD  AND  NEW  TANGENTS" *F 12 .5) 

540  F0RMAT!9H  CURVCF  ■tE14.Ai9H  CURVXF  x,Ei4.4) 

570  FORMAT!  14H  USING  STEPSIZE=>.F12 ,5) 

580  F0RNAT!12H  PREDICTED  X/lXi5F12.5) 

590  FGRNAT!1H  *t2d4H  STEP  REDUCTIONS) 

400  F0RMAT!7H  0MEGA=.F12.5»7H  THETA** FI 2. 5»9H  NEW  RAD«*F12.5) 

410  F0RNAT!31H  TANGNT  ANTICIPATES  I.INTT  POINT) 


PTCM124J 
PTCN1247 
PTCNI243 
PTi:N1244 
PTCN1245 
PTt:N12A4 
PTCN12A7 
PTCN1248 
PTCN12A9 
PTCN1770 
PTCN1271 
PTCN1272 
PTCN1273 
PTCN1274 
PTCN12;5 
PTCH1274 
PTCN1277 
PTCN1278 
PTCN1279 
PTCN1280 
PTCN1281 
PTCN1282 
PTCN1283 
PTCN1284 
PTCN1285 
P TON 1284 
PTCN1287 
PTCN1288 
PTCM1289 
PTCN1290 
PTCN1291 
PTCN1292 
PTCM1293 
PTCN1294 
PTCN1295 
PTCN1294 
PTCMJ297 
PTCN1298 
PTCN1299 
PTCNt300 
PTCM1301 
PTCN1302 
PTCNJ303 
PTCN1304 
PTCNI305 
PTCN130A 
PTCN1307 
PTCN1.308 
PTCN1309 
PTCN1310 
PTCN1311 
PTCN1312 
PTCN13J3 
PTCN  13.14 
PTCNniS 
PTCH131A 
PTCN1317 
PTCN1318 
PTCNI319 
PTCN 1320 


C  PTCN1321 

C  PTCN1323 

END  PTCH1324 
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SUBROUTINE  CORECT  <  NVAR » X » IHOl.D  *  WORK  •  IF.RR .  IMOU*  FPRYM ,  IPVT ,  CRCTOOOI 

1  ABSERR»REI.ERR»XSTEP>NE0N»FNRH)  rRCTOOOr 

C  CRCTOOOS 

•  CRCTOOOS 

SUBROUTINE  CORECT  PERFORMS  THE  CORRECTOR  ITERATIONS  ON  A  STARTING 
POINT.  THE  CORRECTION  METHOD  IS  EITHER  FULl.  (IM0D=0) 

OR  MODIFIED  (IMOD=l)  NEUTON'S  METHOD.  FOR  MODIFIED  NEWTON'S 
METHOD*  THE  JACOBIAN  IS  TO  BE  EUALUATED  ONLY  AT  THE  STARTING  POINT. 

IF  B  IS  THE  VAI.UE  OF  X( IHOl.D)  FOR  THE  INPUT  STARTING  POINT* 

THEN  THE  AUGMENTING  EQUATION  IS  X(IHni.D)=B* 

THAT  IS*  THE  IHOLD-TH  COMPONENT  OF  X  IS  TO  BE  HELD  FIXED. 

THE  AUGMENTED  SYSTEM  TO  BE  SOI.VED  IS  THEN  nFA(X*IHOl.D)*nELTA=FA(X) 


INPUT 

X  =  THE  STARTING  POINT  FOR  THE  CORRECTOR  ITERATION. 

.THOin  s  COMPONENT  OF  X  THAT  WILL  NOT  BE  CHANGED  DURING  ITERATION 
IMOD  *  FLAG  FOR  TYPE  OF  NEWTON'S  METHOD  TO  BE  USED. 

WHEN  IH0D»O*  .JACOBIAN  IS  TO  BE  EVALUATED  AT  EVERY 
CORRECTOR  ITERATE.  KNMAX  IS  SET  TO  10 
IF  IH00>1*  THE  .JACOBIAN  IS  ONLY  EVALUATED  AT  THE  STARTING 
POINT*  AND  KNMAX  IS  SET  TO  20. 


OUTPUT 
X 

WORK 
lERR 


KN 


=  SOLUTION  VECTOR  ON  A  SUCCESSFUL  CALI.  TO  CORECT. 
a  THE  RESIDUAL  F<X)*  AFTER  A  SUCCESSFUL  CALL  TO  CORECT. 
a  THE  RETURN  FLAG  WITH  THE  FO).LOWING  VALUES 
-2  MAXIMUM  NUMBER  OF  CORRECTOR  ITERATIONS  HERE  TAKEN. 
-1  ERROR  RETURN  FROM  SOl.VE  CALLED  BY  CORECT. 

0  SUCCESSFUL  CORRECTION.  VECTOR  X  RETURNED  SATISIFES 
ABS<F<X)).LE.ABSERR 

1  CORRECTOR  STEP  WAS  UNACCEPTABLE*  CORRECTION  FAILED, 
a  THE  NUMBER  OF  CORRECTOR  ITERATIONS  TAKEN  ON  THIS  CAL.L 


THIS  iUIBROUTINE  IS  CALLED  BY 
PITCON 
AND  CALLS 
SOLVE 

FORTRAN  ABS 
LINPAK  ISAMAX 
LINPAK  SAXPY 
USER  FCTN 


crctoaoa 

CRCT0007 

CRCTOOOS 

CRCT0009 

CRCTOOlO 

CRCTOOll 

CRCrOOt’ 

CRCT0013 

r.RCTOOH 

CRCT0015 

CRCTOOlii 

CRCT0017 

CRCT0013 

CRCT0019 

CRCT0020 

CRCT0021 

CRCT0022 

CRCT002S 

CRCT0024 

CRCT002S 

CRCT0026 

CRCT0027 

CRCT0028 

CRCT0029 

CRCTOO.TO 

CRCT0031 

CRCT0032 

CRCT0033 

CRCT0034 

CRCT0035 

CRCT003O 

CRCT0037 

CRCT0038 

CRCT0039 

CRCTft040 

CRCT0041 

CRCT0042 


CRCT0043 
CRCT0045 

REAL  X(NVAR) *UORK(NVAR) *FP«YN(NVAR*NVAR> 

INTEGER  IPVTINVAR) 

COMMON  /COUNTl/  ICRSL* ITN8L*NSTCR*NCNCR*MTRCR*NL.MCR*NLMRT 
COMMON  /COUNT?/  IFEVAL  *IPEVAL.*1S0I.VE*NRED*NRDSUM*KN*KNSUM 
COMMON  /OUTPUT/  IHRITE 
COMMON  /TOL/  EPHACH*EPSATE*EPSQRT 


INITIALIZE 

KNaO 

KNMAXalO 

IF<IN0D.Ca.])KNMAXa2Q 

IERR-0 

FMPaZ.O 

ICALLal 

XSTEP-O.O 

CALL.  FCTN(NVAR*X*HORK) 
IFEVALalFEVALtl 
IMAXalSAMAX (NEON  *  WORK  *  1 ) 
FNRNaABS(UORKCTMAX)) 

WORK (NVAR) -0.0 

STRICTER  ABSERR  TEST  ON  STARTING  POINT 
IF  (FNRM.LE.O.SkABSERR)  GO  TO  40 


CRCT0046 

CRCT0047 

CRCT0048 

CRCT0049 

CRCTOOSO 

CRCT005L 

CRCT00S2 

CRCT0053 

CRCT0054 

CRCT00S5 

CRCT0056 

CRCT00S7 

CRCT0058 

CRCT0059 

CRCT0040 

CRCT0041 

CRCT0042 

CRCT0043 

CRCTOOM 

CRCT00A5 

CRCT00A4 

CRCT00A7 

CRCTOOAS 

CRCT00A9 

rRCT0070 


ooo  nnn  ooo  or>o  ooo  ooo  ooo 


54 


CRCT007t 

ITFRATION  LOOP  CRCT0072 

CRCT0071 

00  20  I»l»  KNMAX  CRCT0074 

KN=I  CRCT007S 

CALL  S01.WF  ( WAR .  X  t  UORK  >  IHOLO » DFTA » IFXP .  lERR .  IC ALl. » IMOO  •  FPRYH » CRCTOO  76 


I  IPVT)  CRCTOO?? 

ICRSl-sICRSL+l  CRnT0078 

IF(IN0n.FQ.l)ICAI.L>0  CRCT0079 

rSOLVEarSOLVF+l  CRCT0080 

IF  (lERR.NE.O)  00  TO  50  CRCTOORl 

FNRHL«FNRH  CRCT0082 

XSTEPL-XSTEP  CRCT0083 

CALL  SAXPY<WAR>-t.0>ynRKfl>X>l)  CRCT00A4 

IMAX-I8AI1AX(NVAR>W)RK>1>  CRCTOORS 

XRTFP>ABI8(U0RK(IHAX))  CRCT00R6 

IHAX>ISAHAX(WAR>Xfl)  CRCT00R7 

XNRH>ABS<  X  < INAX ) )  CRCTOORS 

CALL  FrTN(WAR>XfyORK)  CRCT00R9 

IFEVAL«rFEVAl.41  CRCT0090 

IMAX*ISAMAX(MF0M<U0RKt 1 )  CRCT009\ 

FNRN-ARS(U0RK(IMAX))  CRCT0092 

UORK(WAR)«0.0  CRCT0093 

CRCT0094 

ACCEPTANCE  TEST  CRCT0095 

CRCT0096 

IF  (FNRH.I  E.EPSATE>  GO  TO  60  CRCT0097 

IF  (FNRN.GT.ABSERR)  GO  TO  10  CRCT0098 

IF  (XSTEP.LE.(ABSFRR4RELERR«XNRH))  GO  TO  60  CRCT0099 

CRCTOlOO 

REJECTION  TEST  .  CRCTOlOl 

CRCT0102 

10  IF  (KN.GT.l.AND.XSTEP.GT.IFNPtXSTFPl))  GO  TO  30  CRCT0J03 

IF  <FNRH.GT.<FHP*FNR«L))  GO  TO  30  CRCT0104 

20  FHPal.05  CRCT0105 

60  TO  40  CRCT0.106 

CRCT0107 

(INSOCCESSFWL  STEP  CRCT0108 

CRCT0109 

30  IERR*-3  CRCTOltO 

IF<IMRITE.EQ,2)«RITE(6fl20)  CRCTOHl 

GO  TO  70  CRCT0n2 

CRCT0113 

NAXIWJN  NUNBER  OF  CORRECTOR  STEPS  REACHED  r  .  r,T0X14 

LRCTOI 15 

40  rERR*-2  CRCT0116 

IF<I«RITF,En.2)WRITF(6*110)  CRCT0117 

GO  TO  70  CRCT0118 

CRCT0119 

ERROR  RETURN  IN  SOLVE  CRCT0120 

CRCT0121 

50  IERR—1  CRCT0122 

IF(IMRITE.En.2)«RITE<6flOO)  CRCT0123 

GO  TO  70  CRCT0124 

CRr.T0125 

SUCCESSFUL  STEP  CRCT0126 

CRCT0127 

60  IERR-0  CRCT0128 

70  KNSUN-KNSUN+KN  CRCT0129 

IF(IHRITE»EQ.2)MRITE<6»80)  KN»XSTEP  CRCT0130 

IF(IURITE.Ea.2)URITE<6»90)IH0LD  CRCTOISI 

RETURN  CRCT0132 

80  F0RHAT<13H  CORECT  TOOK  tI2»21H  STEPS.  LAST  ONE  UAS  .E12.5)  CRCT0133 

90  F0RHAT<X4H  CORECT  IHOLD-.H)  CRCT0134 

100  FORMAT ( 31 HOSOLVE  FAILED.  CAI.LED  BY  CORECT)  CRCT0135 

no  FORMAT <25H0T00  MANY  CORRECTOR  STEPS)  CRCTOISA 

120  FORMAT (24HOCORRECTOR  STEP  REJECTED)  CRCT0137 

C  CRCT0t38 

C  CRCT0140 

END  CRCTOIAl 


i 


55 


SIJBROJJTIME  TAN6NT(NVAIl>X»rP,TANf.tRET>ICAl.l.»FPRY/1»IPVT.<<EQN.riETA.  TNr-NOOOl 
1  lEXP)  TNG><0002 

C  TN0N0OO3 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


TNGNOOOS 
TNGNOOOA 
TNGN0007 
TNGN0009 
TNGN0009 
TNGMOOtO 
TMGNOOl  1 

THE  MVAR  BY  MVAR  MATRIX  WHOSE  FIRST  NVAR-.t  ROUSTHGNOO.t2 
DERIVATIVE  OF  FX  EVALUATED  AT  X>  AND  WHOSE  LAST  TNGNOOl.T 


SUBROUTINE  TANGNT  COMPUTES  THE  UNIT  TANGENT  VECTOR  TO  THE  SOLUTION 
CURVE  OF  THE  UNDERDETERHINED  NONLINEAR  SYSTEM  FX  =  0,  THE 
TANGENT  VECTOR  TAN  IS  THE  SOI.UTION  OF  THE  LINEAR  SYSTEM 

DFA<X.IPl.)«TAN  »  E<NVAR) 

WHERE  DFA(X>IPt.)  IS 
ARE  DFX/DX  (X)>  THE 
ROW  IS  <E<IPl.)>  TRANSPOSE*  THE  NVAR  COMPONENT  EUCLIDEAN  COORDINATE 
^CTOR  WITH  I  IN  THE  IPl.-TH  POSITION  AND  ZEROS  ELSEWHERE.  E(NVAR) 


THE  NVAR  COMPONENT  EUCLIDEAN  COORDINATE 
LAST  COMPONENT. 

THE  TANGENT  VECTOR  IS  THEN  NORMALIZED 


VECTOR  WITH  ONE  IN  THE 


AND  ITS  SIGN  ADJUSTED. 


INPUT 
NVAR  > 
X 

IP  * 

OUTPUT 

TAN 

BETA  a 
lEXP  * 
IP  s 


THE  NUMBER  OF  VARIABLES 

THE  CURRENT  CONTINUATION  POINT 

CONTINUATION  COMPONENT  SET  ON  LAST 


STEP 


THE  UNIT  TANGENT  VECTOR  IN  CONTINUATION  DIRECTION  AT  X 
BINARY  MANTISSA  OF  DETERMINANT  OF  JACOBIAN  DFA<X>IPL) 
BINARY  EXPONENT  OF  THE  DETERMINANT  OF  .JACOBIAN  DFA(X*IPL) 
LOCATION  OF  LARGEST  COMPONENT  OF  TANGENT  VECTOR  TAN 
CANDIDATE  FOR  NEW  CONTINUATION  COMPONENT 


THIS  SUBROUTINE 
PITCON 
AND  CALLS 
solve 

LINPAK  ISANAX 
LINPAK  SNRM2 
LINPAK  SSCAL 


IS  CALL.ED  BY 


C 

C 

C 


REAL  XLNVAR) *TAN(NVAR) >FPRYN< NVAR* NVAR) 

INTEGER  IPVT(NVAR) 

COMMtW  /COUNTI/  ICRSL.*ITNSl.*NSTCR*NCNCR*NTRCR*NLMCR*m.MRT 
COMMON  /C0UNT2/  IFFVAL *IPEVAL.*ISOI.VE*NRED*NROSUM*KN*KNSUM 
COMMON  /OUTPUT/  IWRITE 

COMPUTE  TANGENT  VECTOR 


10 


C 

C 

C 


DO  10  I*1*NE0N 
TAN(I)«0.0 
TAN(NVAR)>1.0 
IERR>0 

CALL  SOl.VE<NVAR*X>TAN*IPfDETA*IEXP*IERR*ICALL*IMOn*FPRYN* 
1  IPVT) 

ITNSL-ITNSL41 

IS0LVE>IS0LVE4^1 

IF  (lERR.NE.O)  IRET«<] 

IF  LIRET.LT.O)  RETURN 

OBTAIN  EUCLIDEAN  NORN  OF  TANGENT  VECTOR 

IP«ISANAX(NVAR» TAN* 1 ) 

TN0RN«SNRN2(NVAR*TAN* 1 ) 

IF  (TNORM.EO.O.A)  IRET>-? 

IF  LIRET.LT.O)  RETURN 

NORMALIZE  THE  VECTOR 

,SCALER».1..0/TN0RM 

CAi.L  SSCAL(NVAR*SCAL.ER*TAN*n 

RETURN 


TNGNOOH 
ISTNGNOOIS 
TNGN0016 
TNGN0017 
TNGN0018 
TNGN0019 
TNGNOOZO 
TN6N0021 
TNGN0022 
TNGN0023 
TMGN0024 
TN0N0025 
TNGN002& 
TM6N0027 
TNGN0028 
TNRN0029 
TNGNOO.TO 
TNGN0031 
TNGN0032 
TNGN00.T3 
TNGN0034 
TN6N003.S 
TNGN0036 
TNGN00.T7 
TNGN003B 
fMGN0039 


TNGN0041 
TNGN0042 
TNGN004.T 
TNGN0044 
TNGN004S 
TNGN0044 
TMGN0047 
TMGN0048 
TNGN0049 
TNGNOOBO 
TNGNOOSl 
TNGNOOS2 
TN6N0053 
TN6N0aS4 
TMGN00S5 
TNGN00S6 
TMGM0057 
TNGN0Q58 
TN6N00S9 
TNGN0040 
TNGNOOAl 
TNGN0042 
TNGN00A3 
TNGN00A4 
TNGN0A6S 
TNGN00«A 
TNGN00A7 
TNTtNOOAS 
TNGN0069 
TNGN0070 
TM6N0071 
TN0N0072 

C  TNGN0073 

C««tt«llXt*t*tXXtXtXX<«X««««t«tSX«tS«*X««X««tX«XX««t«««t|[««t«««««tXtttXXTNGN0074 

C  TMGM007S 

PMf  TNGN0076 


ooooooooooonoooooooo  nr>rjooor>ooooooor>ooonoor»r>r>oonnooooor>or>oooooor>oorjo 
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SUBROIITINK  ROOT<AtFAfBiFIt>C>FC>KOIJNTiIFLAn)  ROOTOOOl 

C  RO0T0O02 

R00TrtOO4 

SUBROUTINE  ROOT  SEEKS  A  ROOT  OF  THE  EQUATION  F<X)=0.0>  ROOTCOOS 

OTUEN  A  STARTING  INTERVAL  (AfC)  ON  UHICH  F  CHANGES  SIGN.  R00T0006 

ON  FIRST  CALL  TO  ROOT.  THE  INTERVAL  AND  FUNCTION  VALUES  R00T0007 

FA  AND  FC  ARE  FED  IN  AND  AN  AFPROXIHATION  B  FOR  THE  ROOT  IS  REKJRNED.ROOTOOOS 


BEFORE  EACH  SUBSEQUENT  CALL.  THE  USER  EVALUATES  FB«F(B).  AND  THE 
PROGRAN  TRIES  TO  RETURN  A  BETTER  APPROXIHATION  P. 

THIS  PROGRAM  IS  BASED  ON  THE  FORTRAN  FUNCTION  ZERO 
GIVEN  IN  THE  BOOK: 

'ALGORITHMS  FOR  MINIMIZATION  WITHOUT  DERIVATIVES' 

BY  RICHARD  P.  BRENT.  PRENTICE  HALL. INC.  1973 

THE  MODIFICATIONS  HERE  DONE  BY  JOHN  BURKARDT. 


ON  input: 


A 

FA 


FB 


FC 


IS  ONE  ENDPOINT  OF  AN  INTERVAL  IN  WHICH  F  CHANGES  SIGN. 


ROOT0OO9 
ROOTOOl 0 
ROOTOOll 
ROOTOOl r 
R00T0Ot3 
ROOTOOl 4 
ROOTOOl 5 
R0OT0016 
ROOTOOl 7 
ROOTOOl 8 
ROOTOOl 9 
ROOT0020 
RO0T0021 


-  THE  VALUE  OF  F<A).  THE  USER  MUST  EVALUATE  F(A)  BEFORE  FJRSTRO0T0022 


CALL  ONLY.  THEREAFTER  THE  PROGRAN  SETS  FA. 

ON  FIRST  CALL.  B  SHOULD  NOT  BE  SET  BY  THE  USER. 

ON  SUBSEQUENT  CALLS.  B  SHOULD  NOT  BE  CHANGED 
FROM  ITS  OUTPUT  VALUE.  THE  CURRENT  APPROXIMANT 
TO  THE  ROOT. 

ON  FIRST  CALL.  FB  SHOULD  NOT  BE  SET  BY  THE  USER. 
THEREAFTER.  THE  USER  SHOULD  EVALUATE  THE  FUNCTION 
AT  THE  OUTPUT  VAL.UE  B.  AND  RETURN  FB»F(B). 

IS  THE  OTHER  ENDPOINT  OF  THE  INTERVAL  IN  UHICH 
F  CHANGES  SIGN,  NOTE  THAT  THE  PROGRAM  WILL  RETURN 
IMMEDIATELY  WITH  AN  ERROR  FLAG  IF  FCtFA.GT.0.0. 


R00T0023 

R00T0024 

R00T0025 

RO0T0O26 

R00T0027 

R00T0028 

ROOT0O29 

RO0T0O30 

R0OT0031 

R00T0032 

R00T0033 


THE  VAI.UE  OF  F<C).  THE  USER  MUST  EVALUATE  F(C)  BEFORE  FIRSTROOT0034 


KOUNT 


IFLAG 


CALL  ONLY.  THERAFTER  THE  PROGRAM  SETS  FC. 

-  A  COUNTER  FOR  THE  NUMBER  OF  CAl.l.S  TO  ROOT.  KO:JNT 
SHOULD  BE  SET  TO  ZERO  ON  THE  FIRST  CALL  FOR  A  GIVEN 
ROOT  PROBLEM. 

-  AN  ERROR  RETURN  FLAG  WHOSE  INPUT  VALUE  IS  IMMATERIAL. 


ON  RETURN  FROM  A  CALL  TO  ROOT 


A 

FA 

B 

FB 

C 

FC 

KOUNT 

IFLAG 


-  ONE  ENDPOINT  OF  CHANGE  OF  SIGN  INTERVAL. 

-  THE  VMUE  OF  F<A). 

-  CURRENT  ^PROXIMATION  TO  THE  ROOT.  BEFORE  ANOTHER 
CALL  TO  ROOT.  EVALUATE  F(B). 

-  FB  WILL  BE  OVERWRITTEN  BY  THE  USER  BEFORE  ANOTHER 
CAl.i..  ITS  VALUE  ON  RETURN  IS  ONE  OF  FA.  FP  OR  FC. 

-  OTHER  ENDPOINT  OF  CHANGE  IN  SIGN  INTERVAL. 

-  THE  VALUE  OF  F<C). 

-  CURRENT  NUMBER  OF  CALLS  TO  ROOT. 

-  PROGRAM  RETURN  FLAG! 

IFLA6»-2  MEANS  THAT  ON  FIRST  CALL.  FAXFC.GT.0.0. 

THIS  IS  AN  ERROR  RETURN.  SINCE  A  BRACKETING 
INTERVAL  SHOULD  BE  SUPPLIED  ON  FIRST  CALL. 
IFLAG»-1  MEANS  THAT  THE  CURRENT  BRACKETING  INTERVAL 
WHOSE  ENDPOINTS  ARE  STORED  IN  A  AND  C 
IS  SO  SMALL  (LESS  THAN  4*EPMACH»ABS(B>(-EPMACH) 
THAT  B  SHOULD  BE  ACCEPTED  AS  THE  ROOT. 

THE  FUNCTION  VALUE  F(B)  IS  STORED  IN  FB. 

IFLAG-  0  MEANS  THAT  Tl^  INPUT  VALUE  FB  IS  EXACTLY 

ZERO.  AND  B  SHOULD  PE  ACCEPTED  AS  THE  ROOT. 


R00T003S 

ROOTO036 

ROOT0O37 

ROOT0038 

RO0T0O39 

R00T0040 

ROOTOOAl 

R00T0042 

R0OT0O43 

R00T0044 

R00T0045 

ROOT0046 

R00T0047 

R00T0048 

ROOT0O49 

ROOTOOSO 

R00T0051 

R00T00S2 

R00T0053 

R00T0054 

ROOT0O5S 

ROOTOO!56 

ROOTOOS7 

ROOTOOS8 

R00T00S9 

ROOT0O6O 

ROOTOOAl 

R00T0062 

ROOTOOA3 

R00T00<S4 

R00T0065 

ROOTOOA4 


IFLAG.GT.O  MEANS  THAT  THE  CURRENT  APPROXIMATION  TO 
THE  ROOT  IS  CONTAINED  IN  B.  IF  A  BETTER 
APPROXIMATION  IS  DESIRED.  SET  FB-F<B) 

AND  CALL  ROOT  AGAIN.  THE  VALUE  OF  IFLAG  INDICATES  ROOTOOA7 
THE  METHOD  THAT  WAS  USED  TO  PRODUCE  B.  ROOTOOAS 

ROOTOOA9 

IFLAG-  1  BISECTION  WAS  USED.  R00TDO70 


ooooo  ofiooo  ooo  ooo  oc»oor>  oooooonoooooooooooooonooooooo 
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IFIA6>  2  LINEAR  INTERPni.ATinN  (SECANT  METHOD). 
IFLAO«  S  INVERSE  QUADRATIC  INTERPOLATION. 


LOCAl.  VARIABLES  INa.()DE: 


ROOT0O7t 

R00T0072 

R00T0073 

R00T0O74 

ROOT0O7S 

R00T0076 

R00T0077 

R0OT0O78 

R00T0079 

R00T0080 

ROOTOOB) 


EPNACH*  SMALLEST  POSITIVE  NimHER  SUCH  THAT  1 .04EPNACH.BT.1 .0 
.5«DETA««(]-TAU)  FOR  ROUNDED*  TAU-DIGIT  ARITHMETIC 
BASE  BETA.  TWICE  THAT  VA1.(JE  FOR  TRUNCATED  ARITHMETIC. 

THIS  IS  THE  RELATIVE  MACHINE  PRECISION. 

HALFBC-  SI(»IED  HALFUIDTH  OF  INTERVAL.  DURING  SEGMENT  3*  THE 

CHANGE  OF  SIGN  INTERVAL  IS  <8*0  OR  (C>B>.  THE  MIDPOINT 
OF  THAT  INTERVAL  IS  XNID>Bf HALFBC*  REGARDLESS  OF  ORIENTATION. R00T00B2 
SDELl  -  SIZE  OF  CHANGE  IN  SIGN  INTERVAL.  R00T00B3 

SDELZ  -  PREVIOUS  UAI.UE  OF  SDELl.  R00T0084 

SDEL3  -  PREVIOUS  VALUE  OF  SDEL2.  R00T0085 

SnEL4  -  PREVIOUS  VAI.UE  OF  SDEL3.  R00TOO86 

STEP  -  THE  NEW  ROOT  IS  COMPUTED  AS  A  CORRECTION  TO  B  OF  THE  RO0T0087 

FORM  B<NEU)=B<OLD)+STEP.  R00T0088 

TOI.ER  -  A  NUMBER  WE  ACCEPT  AS  'SMAl.L'  WHEN  EXAMINING  INTERVAL  R00T00B9 

SIZE  OR  STEP  SIZE.  T0LER»2.0»EPMACH*ABS(B)  f  EPMACH  IS  R00T0090 

A  MINIMUM  BELOW  WHICH  WE  HILL  NOT  ALLOW  SUCH  VALUES  TO  FAl.L.  R00T009t 
THIS  SUBROUTINE  IS  CALLED  BY  ROOT0O92 

PITCON  R00T0093 

AND  CALLS  R0OTOO94 

FORTRAN  ABS  R00TO095 

FORTRAN  SIGN  ROOT0O9A 

R0OTO097 

R00T0099 

REAL  A*BiCtFA*FBfFCfSTEP*T(».ER*P*G*RfS 
COMMON  /TOL/  EPMACH >EPSATE»EPSQRT 


SECWIENT  i:  FIRST  CALL  HANDLED  SPCCIA1.LY.  DO  BOOKKEEPING. 

SET  CERTAIN  VALUES  ONLY  FOR  INITIM.  CAl.L  WITH  K0UNT«0 

TF  (ICOUNT.GT.O)  60  TO  10 
IF  <FA.GT.0.0.AND.FC.6T.(U0)  BO  TO  110 
IF  (FA.LT.O.O.ANO.FC.LT.O.O)  60  TO  110 
K0UMT»1 

SDEL1«2.0«ABS(C-A) 

S0EL2«2.0»S0EL1 

S0CL3>2.0«SnEL2 

FB-FC 
60  TO  20 

ON  EVERY  CALL*  INCREMENT  COUNTER 
10  KOUNT«KOUNT>] 

RETURN  IF  HIT  MACHINE  ZERO  FOR  F<B) 

IF(FB.EQ.O.O)  GO  TO  90 

SEGMENT  2:  REARRANGE  POINTS  AND  FUNCTION  VALUES  IF 
NECESSARY  SO  THAT  FBBFC.LT.O.O*  AND  SO  THAT 
ABS(FB) .LT.ABS<FC) 

IF((FB.LE.O.O).AND.(FC.GT.O.O>>  SO  TO  30 
IF((FB.6T.O.O).AND.(FC.LE.O.O>)  GO  TO  30 

FB  AND  FC  ARE  OF  SANE  SIGN. 

(ROOT  CHANGED  SIGN) 

OVERWRITE  C  WITH  VALUE  OF  A 

20  C>A 
FC«FA 


ROOTOIOQ 

ROOTOlOl 

R00T0102 

R00T0103 

R00T0104 

R00T01O5 

R00T0104 

R0OT0107 

R00T0108 

R00T0J09 

ROOTOnO 

ROOTOin 

R00T01.t2 

R00T0113 

R00T0114 

R00T0115 

R00T0116 

R00T0117 

ROOTOnS 

R00T0U9 

R00T0120 

ROOT0121 

R0OT0122 

ROOT0123 

R00T0124 

R0OTO12S 

ROOT0126 

ROOTOJ27 

R00T0128 

ROOT0129 

R00T0I3A 

R00T0131 

R00T0132 

R0OT0133 

R00T0134 

R00T013S 

R00T0136 

R00T0137 

ROOT0138 


ooonooo  or>o  ooo  onoo  oor>o  oooooori  c»r»oo  oooo 


IF  necessaryf  sft  A:«Rf  b:*Cf  c:-b 

TO  ENSURE  THAT  ABS(FB)  .l.E.ABStFC) 

M  TF(ABS(FC).RE.ABS<FB))  60  TO  40 
A*B 
B*C 
C=A 
FA«FB 
FB=FC 
FCaFA 

SEGMENT  3:  CHECK  FOR  ACCEPTANCE  BECAUSE  OF  SHALL  INTF.ROAl. 
CURRENT  CHANGE  IN  SIGN  INTERVAL  IS  (C>B)  OR  (BfC). 

40  T0l.ER-2.0«EPHACH«ABS(B>4-EPHACH 
HAI.FBC=0.5*<C-B) 

SnEL4*SDEL3 

SnEL3»SDEL? 

SDE1.2«SDEL1 

SDEL1*ABS(C-B) 

IF(ABS(HALFBC).t.E. TOLER)  GO  TO  100 

SEGMENT  4*.  COMPUTE  NEW  APPROXTMANT  TO  ROOT  OF  THE  FORM 
B<NEH)>B<QLIl)FSTEP. 

METHODS  AVAILABLE  ARE  LINEAR  INTERPOLATION 
INVERSE  OIJADRATIC  INTERPOLATION 
AND  BISECTION. 

IF<ABS(FB).GE.ABS(FA>)GO  TO  70 
IFIA.NE.C)  GO  TO  SO 

ATTEMPT  LINEAR  IMTERP01.ATI0M  IF  ONL  Y  TWO  POINTS  AVAIl.ABl.E 
COMPUTE  P  AND  Q  FOR  APPROXIMATION  B(NEM)=B(0LD)+P/0 


IFl.A6»2 

S-FB/FA 

P>2.0«HALFBC»S 
Q=1,0-S 
GO  TO  60 

ATTEMPT  INVERSE  QUADRATIC  INTERPOLATION  IF  THREE  POINTS  AVAILABLE 
COMPUTE  P  AMO  0  FOR  APPROXIMATION  B(NEH)»B(0LD)+P/0 

SO  IFLA6«3 
S»FB/FA 
Q»FA/FC 
R-FB/FC 

PsS«  ( 2 .  OBHAI.FBCBOB  ( 0-R )  -  <  B’-A  )  « ( R- 1 . 0  )  > 
0»<0-l,0)*(R-1.0)«<S-1.0) 

CORRECT  THE  SIGNS  OF  P  AND  G 

40  IF<P,6T.0.0)a«-Q 
P»ABS<P) 

IF  P/Q  IS  TOO  LARGE*  GO  BACK  TO  BISECTION 

IF(8,0«S0EL1,0T,SDEL4)  GO  TO  70 
IF  (P.8E.].S«ABS(HALFBC«0)-ABS<T0I.ER«0>)  GO  TO  70 
STEP«P/Q 
60  TO  SO 

PERFORM  bisection: 

IF  ABS<FB).6E.AB8<FA) 

OR  INTERPOI.ATION  IS  UNSAFE  (P/0  IS  LARGE) 

OR  IF  THREE  CONSECUTIVE  STEPS  HIVE  NOT  DECREASED 
THE  SIZE  OF  THE  INTERVAL.  BY  A  FACTOR  OF  8,0 

70  IFl.AG-J 
STEP«HAL.FBC 
60  TO  SO 


R00T0139 

R0OT014O 

RO0T0141 

ROOT0142 

R00T0143 

R00T0144 

R00T014S 

R00T0144 

R00T0J47 

R00T0148 

R00T0149 

ROOTOISO 

R00T0151 

R00T0152 

ROOT0153 

R00T0154 

ROOTOISS 

R00T0156 

ROOT0157 

ROOTOISB 

R00T0159 

R00T0140 

R00T0141 

ROOTO)42 

R00T01A3 

R00T0144 

R00T0145 

R00T0144 

R00T0147 

ROOTOIAS 

R0OT01A9 

R00T0170 

ROOT0171 

R0OT0172 

ROOT0173 

R00T0174 

R00T017S 

R00T0176 

ROOTOJ77 

R00T017B 

R0OT0179 

ROOT0180 

R0OT0181 

ROOTOt82 

R00T0183 

ROOT0184 

ROOTOISS 

R00T0184 

R00T0187 

ROOT0188 

R00T0189 

R00T0190 

R00T0191 

R00T0192 

R00T0193 

R00T0194 

R00T0195 

R00T0194 

R00T0197 

R00T0198 

R00T0199 

R00T0200 

R00T0201 

R00T0202 

ROOT02O3 

R00T0204 

R00T0205 

R00T0204 

R00T0207 

ROOT020R 

ROOTOrOR 

ROOTO?tO 


ooo  ooooo  noooo  or>r»on 
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SCGUFMT  S:  VALUE  OF  STEP  HAS  BEFN  COHPUTFD. 
UPDATE  inforhatioh:  a:«b>  fas«fb»  b:»b+step. 
CHANGE  TN  SIGN  INTERVAl.  IS  NON  (A»C)  OR  <C>A). 

SO  A>B 
FA«FB 

IF<ABS(STEP) .LE. TOLER)  STEP»SIGN(T01.ER»HALFBr) 

B-B+STEP 

RETURN 

SPECIAL  RETURNS 

INPUT  POINT  B  IS  EXACT  ROOT 

90  IF1.AG«0 


ROOTO^n 
R00T02J  2 
R00T02.t3 
R00T0214 
R0OT021S 
R00T0216 
ROOT0217 
R00T0218 
R0OT0219 
R0OTO220 
R0OT0221 
R00T0222 


R00T0223 


R0OT0224 


R00T0225 


R00T0224 


RETURN 

CHANGE  IN  SION  INTERVM.  IS  OF  SIZE  LESS  THAN  4«EPHACHtABS(B)iEPNACH 
INTERVAL  RETURNED  AS  (B*C)  OR  (C*B). 

ACCEPT  B  AS  ROOT  UITH  RESIDUAL  F<B)  STORED  IN  FB» 

100  in.Ae«-i 
A«B 
FA-FB 
RETURN 


R00T0227 

R0flT0228 

R00T0229 

R00T0230 

R00T02S1 

R00T0232 

R0OT023S 

R0(IT0234 

RO0TO235 

R00T023A 

R00T0237 


CHANGE  OF  SIGN  CONDITION  VIOLATED  R0OT0238 

R00T0239 

110  IFLA0»-2  R00T0240 

KOUNT>0  R00T0241 

RETURN  ROOT0242 

C  R00T0243 

C  R00T0245 

END  R00T0244 
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SUmnUTINE  SDLVE(NVARiX>YfIP>DFTA.IEXP>IERRiJCAl.L>IHODfFPRYN> 
1  IPV1) 


SLVEOOOl 
SLVE0002 

C  Sl.yE0003 

SLVEOOOS 

IS  CALLED  BY 


THIS  SUBROiJTIHE 
CORECT 
TANGNT 
AND  CM.LS 
FORTRAN  ABS 
LJNPAK  S6EFA 
LINPAK  S6FSL 
USER  FPRIHE 


THIS  SUBROUTINE 
WHERE  DFAIXfIP) 


SOLVES  THE  LINEAR  SYSTEM  DFA(X»IP)*yOUT  =  YIN 
IS  THE  (NVAR)X(NVAR)  MATRIX  WHOSE  FIRST  NUAR  - 
ROUS  ARE  THE  JACOBIAN  COMPUTED  BY  FPRIME>  AND  WHOSE  LAST 
ROW  IS  ALL  0  EXCEPT  FOR  A  I  IN  T}£  IP-TH  COMPONENT. 


1 


YIN  IS  THE  NVAR  COMPONENT  VECTOR  Y  ON  INPUT*  AND  THE  SOLUTION 
VECTOR  TOUT  IS  RETURNED  IN  Y  ON  OUTPUT  AFTER  A  SUCCESSFUL 
SETUP  AND  SOLUTION. 

««NOTE««  SUBROUTINE  SOLVE  USES  FULL  MATRIX  STORAOE  TO  SOI VE  THE 
LINEAR  SYSTEM.  THE  USER  MAY  WISH  TO  REPLACE  THIS  ROUTINE  WITH 
ONE  MORE  SUITED  TO  HIS  PROBLEM. 


BETA 

lEXP 

IMOD 


ICALL 

INFO 

lERR 

NVAR 

X 

Y 

FPRYM 

IPVT 

IP 


BINARY  MANTISSA  OF  THE  DETERMINANT  OF  JACOBIAN  DFA(X*IP) 
BINARY  EXPONENT  OF  THE  DETERMINANT  OF  JACOBIAN  Iff^A(XiIP) 
NEWTON  METHOD  FLAG. 

IMOPaO*  .lACOBTAN  IS  TO  BE  EVA1.UATED  FOR  EVERY  CORRECTOR  STEP 
AND  EVERY  TANGENT  CA1.CULATION 

IMODal*  JACOBIAN  IS  TO  BE  EVALUATED  ONt.Y  FOR  FIRST  CORRECTOR 
STEP*  AND  EVERY  TANGENT  CAl.a«.ATION 
SET  UP  FLAG. 

IF  (ICALL.EO.O.AND.IMOD.NE.O)  DON'T  RE-EVAl.UATED  JACOBIAN 
OUTPUT  FROM  8GEFA.  IF  INFO.NE.O*  SGEFA  FOUND  A  7ERO 
PIVOT  WHEN  ELIMINATING  INFO-TH  VARIABLE. 

RETURN  FLAG*  0  MEANS  SUCCESSFUL  SOLUTION*  I  MEANS  FAILURE 

THE  NfJMBER  OF  VARIABt.ES  IN  THE  NONLINEAR  SYSTEM 

THE  POINT  AT  WHICH  TO  EVALUATE  FPRYH 

THE  RIGHT  HAND  SIDE  ON  INPUT*  THE  SOl.UTION 

ON  OUTPUT 

ARRAY  WHERE  DFA(X*IP)  IS  TO  BE  STORED. 

INTEGER  WORK  SPACE  FOR  PIVOT  ROW  SWITCHES  DEMANDED  BY  SGEFA 
THE  VARIABLE  USED  IN  THE  AUGMENTING  EQUATION  THAT  IS  OF  THE 
FORM  X(IP)3B.  HENCE  THE  LAST  ROW  OF  DFA<X*IP)  IS  ALL 
ZERO  EXCEPT  FOR  A  1.0  IN  THE  IP-TH  COLUMN. 


SLVE0006 
SLVE0007 
SLVE0008 
SLVE0009 
SLVEOOlO 
SLVEOOl ] 
SLVE0012 
SLVEOOn 
SLVEOOl 4 
SLVEOOtS 
SLVEOOi6 
SLVEOOl 7 
SLVEOOl 8 
SLVEOOl? 
SLVE0020 
SLVE0021 
SLVE0022 
SLVE0023 
SLVE0024 
SLVE0025 
SLVE002« 
SLVE0027 
SI.VE0028 
SLVE0029 
SI.VE0030 
SLVE0031 
SI.VE0032 
SLVE0033 
SLVE0034 
81.VE0035 
SLVE0036 
SLVE0037 
SLVE0038 
SLVE0039 
SLVE0040 
SLVE0041 
SLVE0042 
SLVE0043 
SLVE0044 
SLVE0045 
SLVF0040 
SLVE0047 
SLVE0048 


SLVE0049 

tBAMBBtBBBBBBABtMtBBtBBBAMBBmBBlElcAFBMFFFBBFtkFtFFFFFFXtXtttFFFkBBSLVEOOSO 

SLVE0051 

REAL  X(NVAR)*Y(NVAR)*FPRYH(NVAR*NVAR) 

INTEGER  IPVTINVAR) 

COMMON  /C0UNT2/  IFEVAL*IPEVAL *ISOI.VE*NRED*NRDSUH*KN*KNSUM 


C 

c 

c 


DEPENDING  ON  VALUES  OF  ICALL  AND  IMOD*  EITHER  SET  UP 

AUGMENTED  JACOBIAN*  DECOMPOSE  INTO  L-U  FACTORS*  AND  GET  DETERMINANT* 

OR  USE  CURRENT  FACTORED  JACOBIAN  WITH  NEW  RIGHT-HAND  SIDE. 

IF  (ICA1.L.EQ.O.AND.INOD.NE.O)  GO  TO  SO 
CAI.L  FPRIHE(NVAR*X*FPRYN) 

IPEVAL«IPEVAL-M 
DO  10  I»1*NVAR 
10  FPRYH(NVAR*I)«0.0 
FPRyM<NVAR»IP)*1.0 

CARRY  OUT  IN  CORE  LU  DECOMPOSITION  OF  NVAR  BY  NVAR  MATRIX 
AND  USE  PIVOT  INFORMATION  TO  COMPUTE  DETERMINANT 


SLVE0G52 

81.VE00S3 

SLVE00S4 

SLVE0055 

SLVE005A 

SLVE0057 

SLVEOOS8 

SLVE00S9 

S1.VE00(S0 

SLVEOOOl 

SLVE0002 

SI.VE0003 

SLVE0004 

SI.VE0005 

SLVEOOOO 

SI.VE00A7 

SLVE0008 


oonn 


SLVEOOA9 

CALL  SGEFA  <  FPRYM >  NVAR » NVAR  n IPVT  ^ INFO ) 

SIVE0070 

OETA«J .0 

SLVE0071 

lEXP^^O 

SLVE0072 

TU0>2.0 

31VE0073 

00  40  I>]rNVAR 

SLVE0074 

IF  <IPVT<I).ME.I)  DETA*-DETA 

SIVE0075 

nETA»FPRYM<I»I)*nETA 

SIVE007A 

IF  (BETA. EQ. 0.0)  60  TO  AO 

SLVE0077 

20 

IF  (ABS(nETA).GF.l.O)  60  TO  30 

SLVE0078 

DETAaDETABTUO 

StVE0079 

IEXP»IEXP-1 

SLVE0080 

60  TO  20 

SLVE0081 

30 

IF  <ABS(DETA).LT.TIH))  60  TO  40 

SLOE0082 

BETA-DETA/THO 

SIVE0083 

IEXP-IEXP+1 

SLVE0084 

60  TO  30 

SLVE0085 

40 

CONTXNliE 

SIVE008A 

IF  (JNFO.NE.O)  60  TO  AO 

SLUE0087 

SLVE0088 

USIND  L'U  FACTORED  MATRIX •  SOLVE  SYSTEM  USING  FORUARD-BACKUARD 

SLUE0089 

aiNTNATION*  AND  OVERWRITE  RIGHT  HAND  SIDE  WITH  SOLUTION 

SLVE0090 

SLVE0091 

SO 

CALL  S6ESL(FPRYHfNVARfNVARtIPVT»YfO) 

SLVE0092 

IERR«0 

SIVE0093 

RETURN 

SLUE0094 

AO 

IERR«1 

SIVE0095 

INFOsO 

SLVE009A 

RETURN 

SIVE0097 

SLWE0098 

C  SLWEOIOO 

END  SLVEOlOl 
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